FIXME: I believe this does not work.

    Returns a subset of size K from items with the maximum pairwise distance

    References:
        stackoverflow.com/questions/12278528/subset-elements-furthest-apart-eachother
        stackoverflow.com/questions/13079563/condensed-distance-matrix-pdist

    Recurance:
        Let X[n,k] be the solution for selecting k elements from first n elements items.
        X[n, k] = max( max( X[m, k - 1] + (sum_{p in prev_solution} dist(o, p)) for o < n and o not in prev solution) ) for m < n.

    Example:
        >>> # DISABLE_DOCTEST
        >>> import scipy.spatial.distance as spdist
        >>> items = [1, 6, 20, 21, 22]

    CommandLine:
        python -m utool.util_alg --exec-maximum_distance_subset

    Example:
        >>> # DISABLE_DOCTEST
        >>> from utool.util_alg import *  # NOQA
        >>> #items = [1, 2, 3, 4, 5, 6, 7]
        >>> items = [1, 6, 20, 21, 22]
        >>> K = 3
        >>> result = maximum_distance_subset(items, K)
        >>> print(result)
        (42.0, array([4, 3, 0]), array([22, 21,  1]))

def maximum_distance_subset(items, K, verbose=False):
    """
    FIXME: I believe this does not work.

    Returns a subset of size K from items with the maximum pairwise distance

    References:
        stackoverflow.com/questions/12278528/subset-elements-furthest-apart-eachother
        stackoverflow.com/questions/13079563/condensed-distance-matrix-pdist

    Recurance:
        Let X[n,k] be the solution for selecting k elements from first n elements items.
        X[n, k] = max( max( X[m, k - 1] + (sum_{p in prev_solution} dist(o, p)) for o < n and o not in prev solution) ) for m < n.

    Example:
        >>> # DISABLE_DOCTEST
        >>> import scipy.spatial.distance as spdist
        >>> items = [1, 6, 20, 21, 22]

    CommandLine:
        python -m utool.util_alg --exec-maximum_distance_subset

    Example:
        >>> # DISABLE_DOCTEST
        >>> from utool.util_alg import *  # NOQA
        >>> #items = [1, 2, 3, 4, 5, 6, 7]
        >>> items = [1, 6, 20, 21, 22]
        >>> K = 3
        >>> result = maximum_distance_subset(items, K)
        >>> print(result)
        (42.0, array([4, 3, 0]), array([22, 21,  1]))
    """
    from utool import util_decor
    if verbose:
        print('maximum_distance_subset len(items)=%r, K=%r' % (len(items), K,))

    points = np.array(items)[:, None]

    if False:
        # alternative definition (not sure if works)
        distmat = spdist.squareform(spdist.pdist(points, lambda x, y: np.abs(x - y)))
        D = np.triu(distmat)
        remaining_idxs = np.arange(len(D))
        for count in range(len(points) - K):
            values = D.sum(axis=1) + D.sum(axis=0)
            remove_idx = values.argmin()  # index with minimum pairwise distance
            remaining_idxs = np.delete(remaining_idxs, remove_idx)
            D = np.delete(np.delete(D, remove_idx, axis=0), remove_idx, axis=1)
        value = D.sum()
        subset_idx = remaining_idxs.tolist()
        value, subset_idx
        subset = points.take(subset_idx)
        #print((value, subset_idx, subset))

    sortx = points.T[0].argsort()[::-1]
    sorted_points = points.take(sortx, axis=0)
    pairwise_distance = spdist.pdist(sorted_points, lambda x, y: np.abs(x - y))
    distmat = (spdist.squareform(pairwise_distance))

    def condensed_idx(i, j):
        if i >= len(sorted_points) or j >= len(sorted_points):
            raise IndexError('i=%r j=%r out of range' % (i, j))
        elif i == j:
            return None
        elif j < i:
            i, j = j, i
        return (i * len(sorted_points) + j) - (i * (i + 1) // 2) - (i) - (1)

    def dist(i, j):
        idx = condensed_idx(i, j)
        return 0 if idx is None else pairwise_distance[idx]

    @util_decor.memoize_nonzero
    def optimal_solution(n, k):
        """
        Givem sorted items sorted_points
        Pick subset_idx of size k from sorted_points[:n] with maximum pairwise distance
        Dynamic programming solution
        """
        "# FIXME BROKEN "
        assert n <= len(sorted_points) and k <= len(sorted_points)
        if k < 2 or n < 2 or n < k:
            # BASE CASE
            value, subset_idx =  0, []
        elif k == 2:
            # BASE CASE
            # when k==2 we choose the maximum pairwise pair
            subdist = np.triu(distmat[0:n, 0:n])
            maxpos = subdist.argmax()
            ix, jx = np.unravel_index(maxpos, subdist.shape)
            value = distmat[ix, jx]
            subset_idx = [ix, jx]
        else:
            # RECURSIVE CASE
            value = 0
            subset_idx = None
            # MAX OVER ALL OTHER NODES (might not need a full on loop here, but this will definitely work)
            for m in range(k - 1, n):
                # Choose which point to add would maximize the distance with the previous best answer.
                prev_value, prev_subset_idx = optimal_solution(m, k - 1)
                for o in range(n):
                    if o in prev_subset_idx:
                        continue
                    add_value = sum([distmat[o, px] for px in prev_subset_idx])
                    cur_value = prev_value + add_value
                    if cur_value > value:
                        value = cur_value
                        subset_idx = prev_subset_idx + [o]
        return value, subset_idx

    value, sorted_subset_idx = optimal_solution(len(points), K)
    subset_idx = sortx.take(sorted_subset_idx)
    subset = points.take(subset_idx)
    #print((value, subset_idx, subset))
    return value, subset_idx, subset