human
Complete this python program:
def griddata(points, values, xi, method='linear', fill_value=np.nan, rescale=False):
assistant
#Here is the Python program you requested.

    """Interpolate unstructured D-dimensional data.
    Parameters
    points : ndarray of floats, shape (n, D)
    Data point coordinates. Can either be an array of
    shape (n, D), or a tuple of `ndim` arrays.
    values : ndarray of float or complex, shape (n,)
    Data values.
    xi : 2-D ndarray of float or tuple of 1-D array, shape (M, D)
    Points at which to interpolate data.
    method : {\linear\, \nearest\, \cubic\}, optional
    Method of interpolation. One of
    ``nearest``
    return the value at the data point closest to
    the point of interpolation.  See `NearestNDInterpolator` for
    more details.
    ``linear``
    tesselate the input point set to n-dimensional
    simplices, and interpolate linearly on each simplex.  See
    `LinearNDInterpolator` for more details.
    ``cubic`` (1-D)
    return the value determined from a cubic
    spline.
    ``cubic`` (2-D)
    return the value determined from a
    piecewise cubic, continuously differentiable (C1), and
    approximately curvature-minimizing polynomial surface. See
    `CloughTocher2DInterpolator` for more details.
    fill_value : float, optional
    Value used to fill in for requested points outside of the
    convex hull of the input points.  If not provided, then the
    default is ``nan``. This option has no effect for the
    \nearest\ method.
    rescale : bool, optional
    Rescale points to unit cube before performing interpolation.
    This is useful if some of the input dimensions have
    incommensurable units and differ by many orders of magnitude.
    .. versionadded:: 0.14.0
    Notes
    .. versionadded:: 0.9
    Examples
    Suppose we want to interpolate the 2-D function
    >>> def func(x, y):
    ...     return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
    on a grid in [0, 1]x[0, 1]
    >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
    but we only know its values at 1000 data points:
    >>> points = np.random.rand(1000, 2)
    >>> values = func(points[:,0], points[:,1])
    This can be done with..."""
    points = _ndim_coords_from_arrays(points)
    if (points.ndim < 2):
        ndim = points.ndim
    else:
        ndim = points.shape[(-1)]
    if ((ndim == 1) and (method in ('nearest', 'linear', 'cubic'))):
        from .interpolate import interp1d
        points = points.ravel()
        if isinstance(xi, tuple):
            if (len(xi) != 1):
                raise ValueError('invalid    number    of    dimensions    in    xi')
            (xi,) = xi
        idx = np.argsort(points)
        points = points[idx]
        values = values[idx]
        if (method == 'nearest'):
            fill_value = 'extrapolate'
        ip = interp1d(points, values, kind=method, axis=0, bounds_error=False, fill_value=fill_value)
        return ip(xi)
    elif (method == 'nearest'):
        ip = NearestNDInterpolator(points, values, rescale=rescale)
        return ip(xi)
    elif (method == 'linear'):
        ip = LinearNDInterpolator(points, values, fill_value=fill_value, rescale=rescale)
        return ip(xi)
    elif ((method == 'cubic') and (ndim == 2)):
        ip = CloughTocher2DInterpolator(points, values, fill_value=fill_value, rescale=rescale)
        return ip(xi)
    else:
        raise ValueError(('Unknown    interpolation    method    %r    for    %d    dimensional    data' % (method, ndim)))