I'm using a Delaunay triangulation to interpolate in some function values evaluated at a set of parameters on a regular 4-dimensional grid. Sometimes, when a parameter value changes by a small amount that takes it to a new simplex, more than one point in the simplex changes. I expect that as I vary one of the parameter continuously, I'd move from simplex to simplex by changing just one point in the simplex at a time (and that's usually the case in my code, too). Instead, consider this script:
import numpy as np
from scipy.spatial import Delaunay
# hideous construction to get the desired 4d grid of points
# with points at [-1, -0.5, 0, 0.5, 1] along each axis
X = np.vstack(list(map(np.ravel, np.meshgrid(*[np.linspace(-1, 1, 5) for i in range(4)])))).T
tri = Delaunay(X)
delta = 1e-10
print(np.sort(tri.vertices[tri.find_simplex([-0.25, -0.05, 0.5+delta, 0.1])]))
print(np.sort(tri.vertices[tri.find_simplex([-0.25, -0.05, 0.5-delta, 0.1])]))
which produces
[192 292 317 318 322]
[167 292 293 313 317]
Note that these two simplices differ by 3 points, where I expect one, and I haven't devised a 2- or 3-D example where more than one vertex would change.
I'm 99% sure this is because my points are on a regular grid but I can't find a detailed answer of why or how to avoid the problem. I know that the triangulation isn't unique but that isn't fundamentally a problem. Various tricks appear to change where I encounter this issue but I haven't yet found a "fix" that prevents the issue from appearing anywhere.
