Interpolation on n-dimensional grid impossible through scipy.interpolate.griddata

Question

I have an array of observations called values which comes from a function defined the following way:

values = np.array([oscillatory(i) for i in range(points.shape[0])])

values has shape (65,1), points is an array at which the oscillatory function is evaluated, it has shape (65,7) , 7 are some features acting as dimensions in a 7-dimensional space (7 is just an arbitrary number).

I am trying to interpolate some arbitrary points defined in this space. I have tried defining those points with:

grid_x = np.random.uniform(0,10, (100,7))

But it didn't work. Apparently because the grid is not correctly defined so I have tried:

grid_x= np.mgrid[-2:2, 5:99 , 4:5, 5:6, 5:7, 6:67, 7:67]

Which doesn't work again. I am calling the interpolation function the following way:

grid_z1 = gdd(points, values, tuple(grid_x))

But I always get a huge error which I have some trouble to understand.

The curious thing is that if I define points and values in a random way the code works:

points = np.random.uniform(0,10, (65,3))
values = np.random.uniform(0,10,(points.shape[0],1))

grid_x= np.mgrid[0:2, 5:9 , 4:5]
grid_z1 = gdd(points, values, tuple(grid_x))

Here I just tried 3 instead of 7 because its faster, but the principle remains the same. If I define it in 7 dimensions the code also works. So my questions are: 1) How could I get my initial code in 7 dimensions running ? 2) Why is the random declaration working compared to the other if arrays have the same shape ?

Any help would be highly appreciated. Thank you very much.

The error I am getting is the following:

    Traceback (most recent call last):
  File "/usr/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 2, in <module>
  File "/usr/local/lib/python3.6/dist-packages/scipy/interpolate/ndgriddata.py", line 222, in griddata
    rescale=rescale)
  File "interpnd.pyx", line 248, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__
  File "qhull.pyx", line 1828, in scipy.spatial.qhull.Delaunay.__init__
  File "qhull.pyx", line 354, in scipy.spatial.qhull._Qhull.__init__
scipy.spatial.qhull.QhullError: QH6154 Qhull precision error: Initial simplex is flat (facet 1 is coplanar with the interior point)
While executing:  | qhull d Qbb Qt Qz Q12 Qc
Options selected for Qhull 2015.2.r 2016/01/18:
  run-id 526315618  delaunay  Qbbound-last  Qtriangulate  Qz-infinity-point
  Q12-no-wide-dup  Qcoplanar-keep  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood  _max-width  4  Error-roundoff 1.2e-14  _one-merge 1.1e-13
  Visible-distance 7.3e-14  U-coplanar-distance 7.3e-14  Width-outside 1.5e-13
  _wide-facet 4.4e-13
precision problems (corrected unless 'Q0' or an error)
      2 flipped facets
      6 degenerate hyperplanes recomputed with gaussian elimination
     12 nearly singular or axis-parallel hyperplanes
      6 zero divisors during back substitute
    126 zero divisors during gaussian elimination
The input to qhull appears to be less than 4 dimensional, or a
computation has overflowed.
Qhull could not construct a clearly convex simplex from points:
- p2(v4):  -1.9     5     4 0.012
- p1(v3):  -1.9     5     4 0.0054
- p65(v2):     0   5.5   4.5     4
- p64(v1):     2     6     5     3
- p0(v0):    -2     5     4 -8.9e-16
The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 1.2e-14.  The center point, facets and distances
to the center point are as follows:
center point  -0.7625    5.309    4.309    1.407
facet p1 p65 p64 p0 distance=    0
facet p2 p65 p64 p0 distance= -8.9e-16
facet p2 p1 p64 p0 distance= -8.9e-16
facet p2 p1 p65 p0 distance= -8.9e-16
facet p2 p1 p65 p64 distance= -8.9e-16
These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.
The min and max coordinates for each dimension are:
  0:        -2         2  difference=    4
  1:         5         6  difference=    1
  2:         4         5  difference=    1
  3:  -8.882e-16         4  difference=    4
If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 1.2e-14.
  - trace execution with 'T3' to see the determinant for each point.
If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.

Answer 1

I think the problem is that you need to make sure you're consistent about the number of dimensions, and the 'domains' of those dimensions. You won't get good results interpolating into places you have not sampled. I think the errors you're getting are related to trying to compute things in these places.

Here's how to make something like the example in the scipy.interpolate.griddata docs work for a 7-dimensional example. I'm using a much simpler function which just sums the 'features' in the points data:

import numpy as np

def func(data):
    return np.sum(data, axis=1)

grid = np.mgrid[0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j]

points = np.random.rand(100, 7)
values = func(points)

Notice that the grid covers the entire range of the coordinates points. That is, since each column in points has values in the range 0 to 1, I should make sure I'm making a grid over those same coordinates.

from scipy.interpolate import griddata
grid_z = griddata(points, values, tuple(grid), method='linear')

Now I have this:

>>> grid_z[2, 2, 2, 2, 2, :, :]
array([[ nan,  nan,  nan,  nan,  nan],
       [ nan, 3.  , 3.25, 3.5 ,  nan],
       [ nan, 3.25, 3.5 , 3.75,  nan],
       [ nan, 3.5 , 3.75, 4.  ,  nan],
       [ nan,  nan,  nan,  nan,  nan]])

Notice that there are a lot of NaNs. If you use nearest as the method, you'll always get a solution, but of course linear interpolation needs two things to interpolate between, so the 'edges' of the hypercube are not valid (and 7-D space has a lot of edge!).