Appropriate numpy/scipy function to interpolate function defined on simplex (non regular grid)

Question

I've got a function defined on a 3-dimensional simplex. That is, the set of points x, y, z, each between 0 and 1, such that x + y + z = 1.0

For example, if I consider 4 points for each x, y and z, then I'll get a a (10, 3) numpy array that looks like this (each row sums up to exactly 1):

points = array([[0.        , 0.        , 1.        ],
       [0.        , 0.33333333, 0.66666667],
       [0.        , 0.66666667, 0.33333333],
       [0.        , 1.        , 0.        ],
       [0.33333333, 0.        , 0.66666667],
       [0.33333333, 0.33333333, 0.33333333],
       [0.33333333, 0.66666667, 0.        ],
       [0.66666667, 0.        , 0.33333333],
       [0.66666667, 0.33333333, 0.        ],
       [1.        , 0.        , 0.        ]])

I add the convenience function that generates a simplex:

def generate_simplex_3dims(n_per_dim):
    xlist = np.linspace(0.0, 1.0, n_per_dim)
    ylist = np.linspace(0.0, 1.0, n_per_dim)
    zlist = np.linspace(0.0, 1.0, n_per_dim)
    return np.array([[x, y, z] for x in xlist for y in ylist for z in zlist
                     if np.allclose(x+y+z, 1.0)])

I'll also have values for those points. As an example, let's generate the values like this:

def approx_this_f(x, y, z):
    return 2*x - y + 5*z

values = np.empty(len(points))
for i, point in enumerate(points):
    values[i] = approx_this_f(point[0], point[1],
                         point[2])

My objective is to get an interpolated_f that I can use to evaluate like interpolated_f([0.3, 0.5, 0.2]) or interpolated_f(0.3, 0.5, 0.2) for arbitrary points within the simplex.

I looked through the documentation, but don't understand what is the appropriate interpolator here, given that my grid points are defined on a simplex and that I want to get an interpolated function back.

I tried scipy.interpolate.griddata and it only worked with method='nearest' and this one returns an array of values, but I need an interpolated function. I saw other functions on scipythat return an interpolated function, but seem to only work with regular grids.

Thanks!

---- Example with griddata in case it helps ------

from scipy.interpolate import griddata
xi = generate_simplex_3dims(n_per_dim=20) #Generates lots of points
interpolated_grid = griddata(points, values, xi,
         method='linear') #this fails
interpolated_grid = griddata(points, values, xi,
         method='nearest') #this works, but returns a grid, not a function

The method=linear threw and error, but, more im

Answer 1

Thanks to the answer by @user6655984 I figured out how to do it (thanks!)

I gave it a bit more thought and I'm pretty sure that (thought I'd be happy to be corrected):

1. The simplex domain means it can't have a regular grid (some parts of the space won't have a value)
2. That I need to delete one dimension for the interpolation to work. Since in a 3D simplex the elements need to add up to one, I'll have [0.4, 0.3, 0.3], but I can't have [0.4, 0.3, 0.5], so in fact the function will only have values on a 2D subspace of 3D!

Let's have the same setup as in the question

def approx_this_f(x, y, z):
    return 2*x - y + 5*z

#Uses the function defined in the question
simplex_points = generate_simplex_3dims(10)
values = np.empty(len(simplex_points))
for i, lambda_0 in enumerate(simplex_points):
    values[i] = approx_this_f(lambda_0[0], lambda_0[1],
                         lambda_0[2])

Because of comment number 2, the following doesn't work:

from scipy.interpolate import LinearNDInterpolator

interpolated = LinearNDInterpolator(simplex_points, values)

and throws this error

QhullError: qhull precision warning: 
The initial hull is narrow (cosine of min. angle is 0.9999999999999999).
Is the input lower dimensional (e.g., on a plane in 3-d)?  Qhull may
produce a wide facet. 

So I need to pass the points with one less dimension (i.e., only columns 1 and 2, not the third one):

interpolated = LinearNDInterpolator(simplex_points[:, 0:2], values)

Now we can evaluate it on other points

#Silly code to make the original function take a matrix
def approx_this_f_vec(array):
    res = np.empty(len(array))
    for row in range(len(array)):
         res[row] = approx_this_f(*array[row])
     return res

points_test = src.generate_simplex_3dims(50) #1275 new points
interpolated_vals = interpolated_f(points_test[:, 0:2])
real_values = approx_this_f_vec(points_test)

print((interpolated_vals - real_values).max())

gives 1.77e-15 which means the interpolation did pretty well!