delaunay_plot.py

# BSD 3-Clause License

# Copyright (c) 2018, John H. Williamson
# All rights reserved.

# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:

# * Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.

# * Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.

# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.

# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

import scipy.spatial
import scipy.interpolate
from sklearn.metrics.pairwise import paired_distances

def delaunay_distance_plot(low_d, high_d, metric='euclidean', resolution=200):
    triangulation = scipy.spatial.Delaunay(low_d)

    # remove all exterior points
    # and make sure we have a -> b -> c -> a
    simplices = triangulation.simplices[np.all(triangulation.simplices > 0, 
                                               axis=1)]
    simplices = np.concatenate([triangulation.simplices,
                                triangulation.simplices[:, 0:1]],
                                axis=1)
    edges = np.stack([simplices[:, :-1], simplices[:, 1:]])

    # slice and rearrange to 2 x n_vertices x n_edges
    edges = edges[:, :, :triangulation.simplices.shape[1]]

    # split into either ends of the edges
    low_v1, low_v2 = low_d[edges[0, :, :].ravel()], low_d[edges[1, :, :].ravel()]
    high_v1, high_v2 = high_d[edges[0, :, :].ravel()], high_d[edges[1, :, :].ravel()]

    low_distances = paired_distances(low_v1, low_v2, metric=metric)
    high_distances = paired_distances(high_v1, high_v2, metric=metric)

    # compute centre points of edges and ratio of distances
    # along that edge
    ctrs = 0.5 * (low_v1 + low_v2)
    ratios = low_distances / high_distances

    # interpolate ratios over low-d space
    mx, my = np.mgrid[np.min(low_d[:, 0]):np.max(low_d[:, 0]):resolution * 1j,
                      np.min(low_d[:, 1]):np.max(low_d[:, 1]):resolution * 1j]

    gridded = scipy.interpolate.griddata(
        ctrs, ratios**0.5, (mx, my), method='linear')
    
       
    plt.pcolormesh(mx, my, gridded, cmap='viridis')
    
# e.g.
# delaunay_distance_plot(embedding, data)