Torus topology interconnection interpolation

torus_interpolation.py

import itertools
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

points = np.array([
    [0, 0],
    [0, 2],
    [2, 0],
    [0.5, 0.5]
])

values = np.array([1, 1, 1, 4])

xi = np.array([
    [0.3, 0.2],
    [0.7, 1.8]
])


def make_toroidal(points, values):
    shifts = [
        np.amax(points[:, i]) - np.amin(points[:, i])
        for i in xrange(points.shape[1])
    ]
    directions = np.array([0, 1, -1])
    all_shifts = itertools.product(*[
        directions * shift for shift in shifts
    ])

    toroidal = np.empty((0,) + points.shape[1:])
    for total_shift in all_shifts:
        toroidal = np.append(toroidal, points + np.array(total_shift), axis=0)

    repeated_values = np.tile(values, 3 ** points.shape[1])

    return toroidal, repeated_values


def find_barycentric(p, simplex):
    start = simplex[:-1]
    last = simplex[-1]
    t = np.transpose(np.matrix(start) - last)
    coords = np.linalg.solve(t, p - last)  # Cramer? third thing to think of
    return np.append(coords, 1 - sum(coords))


def do_interpolation(points, values, xi, fill_value=np.nan):
    tri = Delaunay(points)  # first thing to implement

    results = np.zeros(xi.shape[0])

    for i in xrange(xi.shape[0]):
        p = xi[i]
        isimplex = tri.find_simplex(p)  # second thing to implement
        if isimplex == -1:
            results[i] = fill_value
            continue
        simplex_indices = tri.simplices[isimplex]
        vs = np.array(
            [values[x] for x in simplex_indices]
        )
        if np.isnan(vs).any():
            results[i] = fill_value
            continue
        simplex = np.array(
            [points[x] for x in simplex_indices]
        )
        barycentric_coords = find_barycentric(p, simplex)
        results[i] = np.dot(barycentric_coords, vs)

    return results, tri.simplices


def add_plot(id, title, results, simplices):
    plt.subplot(id)
    plt.title(title)
    plt.triplot(points[:, 0], points[:, 1], simplices.copy())
    plt.plot(points[:, 0], points[:, 1], 'o')
    plt.plot(xi[:, 0], xi[:, 1], 'o')
    for point, value in itertools.chain(
            itertools.izip(points, values),
            itertools.izip(xi, results)
            ):
        plt.annotate(float(value), point)


plt.figure(1)
add_plot(211, 'Standard', *do_interpolation(points, values, xi))
points, values = make_toroidal(points, values)
add_plot(212, 'Torus', *do_interpolation(points, values, xi))
plt.show()