Python scipy.interpolate.interp2d() Examples

def _make_interpolation(self):
        """
        creates an interpolation grid in H_0, omega_m and computes quantities in Dd and Ds_Dds
        :return:
        """
        grid2d = np.dstack(np.meshgrid(self._H0_range, self._omega_m_range)).reshape(-1, 2)
        H0_grid = grid2d[:, 0]
        omega_m_grid = grid2d[:, 1]
        Dd_grid = np.zeros_like(H0_grid)
        Ds_Dds_grid = np.zeros_like(H0_grid)
        for i in range(len(H0_grid)):
            Dd, Ds_Dds = cosmo2angular_diameter_distances(H0_grid[i], omega_m_grid[i], self.z_d, self.z_s)
            Dd_grid[i] = Dd
            Ds_Dds_grid[i] = Ds_Dds
        self._f_H0 = interpolate.interp2d(Dd_grid, Ds_Dds_grid, H0_grid, kind='linear', copy=False, bounds_error=False, fill_value=-1)
        print("H0 interpolation done")
        self._f_omega_m = interpolate.interp2d(Dd_grid, Ds_Dds_grid, omega_m_grid, kind='linear', copy=False, bounds_error=False, fill_value=0)
        print("omega_m interpolation done") 

def _interpolate(x1, y1, z1, x2, y2):
    """Helper to interpolate in 1d or 2d

    We interpolate to get the same shape than with other methods.
    """
    if x1.size > 1 and y1.size > 1:
        func = interp2d(x1, y1, z1.T, kind='linear', bounds_error=False)
        z2 = func(x2, y2)
    elif x1.size == 1 and y1.size > 1:
        func = interp1d(y1, z1.ravel(), kind='linear', bounds_error=False)
        z2 = func(y2)
    elif y1.size == 1 and x1.size > 1:
        func = interp1d(x1, z1.ravel(), kind='linear', bounds_error=False)
        z2 = func(x2)
    else:
        raise ValueError("Can't interpolate a scalar.")

    # interp2d is not intuitive and return this shape:
    z2.shape = (y2.size, x2.size)
    return z2 

def _infer_interval_breaks(coord, kind=None):
        """
        Interpolate the bounds from the data in coord

        Parameters
        ----------
        %(CFDecoder.get_plotbounds.parameters.no_ignore_shape)s

        Returns
        -------
        %(CFDecoder.get_plotbounds.returns)s

        Notes
        -----
        this currently only works for rectilinear grids"""
        if coord.ndim == 1:
            return _infer_interval_breaks(coord)
        elif coord.ndim == 2:
            from scipy.interpolate import interp2d
            kind = kind or rcParams['decoder.interp_kind']
            y, x = map(np.arange, coord.shape)
            new_x, new_y = map(_infer_interval_breaks, [x, y])
            coord = np.asarray(coord)
            return interp2d(x, y, coord, kind=kind, copy=False)(new_x, new_y) 

def __call__(self, x, y, dx=0, dy=0, assume_sorted=False):
        """  Wrapper to scipy.interpolate.interp2d which preserves the input ordering.

        Parameters
        ----------
        x: See superclass
        y: See superclass
        dx: See superclass
        dy: See superclass
        assume_sorted: bool, optional
            This is just a place holder to prevent a warning.
            Overwriting this will not do anything

        Returns
        ----------
        array_like: See superclass

        """
        unsorted_idxs = np.argsort(np.argsort(x))
        return super(UnsortedInterp2d, self).__call__(x, y, dx=dx, dy=dy, assume_sorted=False)[unsorted_idxs]


#  Instantiate the default argument parser at runtime 

def selected_interpolation(fin,mx,my,mz,nite=3):
  """Call in a function in a mesh which is two nite times finer, but only
  in those points that are expected to change, in the others interpolate"""
  if nite<1.01: return mx,my,mz 
  f = interpolate.interp2d(mx[:,0], my[0,:], mz, kind='linear') # interpolation
  x = np.linspace(np.min(mx),np.max(mx),len(mx)*2) # twice as big
  y = np.linspace(np.min(my),np.max(my),len(my)*2) # twice as big
  mx2, my2 = np.meshgrid(x, y) # new grid
  mz2 = f(x,y) # interpolate
  mz2 = mz2.transpose()
  dmz = np.gradient(mz2) # gradient
  dmz = dmz[0]*dmz[0] + dmz[1]*dmz[1] # norm of the derivative
  maxdm = np.max(dmz) # maximum derivative
  for i in range(len(mx2)):
    for j in range(len(my2)):
      if dmz[i,j]>0.001*maxdm: # if derivative is large in this point
        mz2[i,j] = fin(mx2[i,j],my2[i,j]) # re-evaluate the function
  mx2 = mx2.transpose()
  my2 = my2.transpose()
  mz2 = mz2.transpose()
  if nite>1:
    return selected_interpolation(fin,mx2,my2,mz2,nite=nite-1)
  return mx2,my2,mz2 # return function 

def interpolator2d(x,y,z):
  """Return a 2D interpolator"""
  x = np.array(x)
  y = np.array(y)
  from scipy.interpolate import interp2d
  ny = (y.max()-y.min())/abs(y[1]-y[0])
  ny = int(round(ny)) + 1
  nx = len(z)/ny
  nx = int(nx)
  ny = int(ny)
  x0 = np.linspace(min(x),max(x),nx)
  y0 = np.linspace(min(y),max(y),ny)
  xx, yy = np.meshgrid(x0, y0)
  Z = np.array(z).reshape(nx,ny) # makes a (Zy,Zx) matrix out of z
  T = Z.T
  f = interp2d(x0, y0, T, kind='linear')
  return f 

def interpol_lin2d(x, y, z, xout, yout):
    """
    Performs a linear interpolation of a 2D matrix/image irregularily sampled on both axes

    :author: Andreas Reigber
    :param x: The x values of the first axis of z
    :type x: 1-D ndarray float
    :param y: The y values of the second axis of z
    :type y: 1-D ndarray float
    :param z: The input matrix
    :type z: 2D ndarray
    :param xout: The values on the first axis where the interpolates are desired
    :type xout: 1D ndarray float
    :param yout: The values on the second axis where the interpolates are desired
    :type yout: 1D ndarray float
    :returns: The interpolated matrix /  image
    """
    if np.iscomplexobj(z):
        f_real = sp.interpolate.interp2d(x, y, z.real, kind='linear')
        f_imag = sp.interpolate.interp2d(x, y, z.imag, kind='linear')
        return f_real(xout, yout) + 1j * f_imag(xout, yout)
    else:
        f = sp.interpolate.interp2d(x, y, z, kind='linear')
        return f(xout, yout) 

def interpol_spline2d(x, y, z, xout, yout):
    """
    Performs a cubic spline interpolation of a 2D matrix/image irregularily sampled on both axes

    :author: Andreas Reigber
    :param x: The x values of the first axis of z
    :type x: 1-D ndarray float
    :param y: The y values of the second axis of z
    :type y: 1-D ndarray float
    :param z: The input matrix
    :type z: 2D ndarray
    :param xout: The values on the first axis where the interpolates are desired
    :type xout: 1D ndarray float
    :param yout: The values on the second axis where the interpolates are desired
    :type yout: 1D ndarray float
    :returns: The interpolated matrix /  image
    """
    if np.iscomplexobj(z):
        f_real = interpolate.interp2d(x, y, z.real, kind='cubic')
        f_imag = interpolate.interp2d(x, y, z.imag, kind='cubic')
        return f_real(xout, yout) + 1j * f_imag(xout, yout)
    else:
        f = interpolate.interp2d(x, y, z, kind='cubic')
        return f(xout, yout) 

def resample_2d(array, sample_pts, query_pts, kind='linear'):
    """Resample 2D array to be sampled along queried points.

    Parameters
    ----------
    array : `numpy.ndarray`
        2D array
    sample_pts : `tuple`
        pair of `numpy.ndarray` objects that contain the x and y sample locations,
        each array should be 1D
    query_pts : `tuple`
        points to interpolate onto, also 1D for each array
    kind : `str`, {'linear', 'cubic', 'quintic'}
        kind / order of spline to use

    Returns
    -------
    `numpy.ndarray`
        array resampled onto query_pts via bivariate spline

    """
    interpf = interpolate.interp2d(*sample_pts, array, kind=kind)
    return interpf(*query_pts) 

def interpolate_xarray(xpoints, ypoints, values, shape, kind='cubic',
                       blocksize=CHUNK_SIZE):
    """Interpolate, generating a dask array."""
    vchunks = range(0, shape[0], blocksize)
    hchunks = range(0, shape[1], blocksize)

    token = tokenize(blocksize, xpoints, ypoints, values, kind, shape)
    name = 'interpolate-' + token

    from scipy.interpolate import interp2d
    interpolator = interp2d(xpoints, ypoints, values, kind=kind)

    dskx = {(name, i, j): (interpolate_slice,
                           slice(vcs, min(vcs + blocksize, shape[0])),
                           slice(hcs, min(hcs + blocksize, shape[1])),
                           interpolator)
            for i, vcs in enumerate(vchunks)
            for j, hcs in enumerate(hchunks)
            }

    res = da.Array(dskx, name, shape=list(shape),
                   chunks=(blocksize, blocksize),
                   dtype=values.dtype)
    return DataArray(res, dims=('y', 'x')) 

def test_interp2d_meshgrid_input_unsorted(self):
        np.random.seed(1234)
        x = linspace(0, 2, 16)
        y = linspace(0, pi, 21)

        z = sin(x[None,:] + y[:,None]/2.)
        ip1 = interp2d(x.copy(), y.copy(), z, kind='cubic')

        np.random.shuffle(x)
        z = sin(x[None,:] + y[:,None]/2.)
        ip2 = interp2d(x.copy(), y.copy(), z, kind='cubic')

        np.random.shuffle(x)
        np.random.shuffle(y)
        z = sin(x[None,:] + y[:,None]/2.)
        ip3 = interp2d(x, y, z, kind='cubic')

        x = linspace(0, 2, 31)
        y = linspace(0, pi, 30)

        assert_equal(ip1(x, y), ip2(x, y))
        assert_equal(ip1(x, y), ip3(x, y)) 

def test_interp2d_bounds(self):
        x = np.linspace(0, 1, 5)
        y = np.linspace(0, 2, 7)
        z = x[:,None]**2 + y[None,:]

        ix = np.linspace(-1, 3, 31)
        iy = np.linspace(-1, 3, 33)

        b = interp2d(x, y, z, bounds_error=True)
        assert_raises(ValueError, b, ix, iy)

        b = interp2d(x, y, z, fill_value=np.nan)
        iz = b(ix, iy)
        mx = (ix < 0) | (ix > 1)
        my = (iy < 0) | (iy > 2)
        assert_(np.isnan(iz[my,:]).all())
        assert_(np.isnan(iz[:,mx]).all())
        assert_(np.isfinite(iz[~my,:][:,~mx]).all()) 

def findpeak3(f,subpixel):
    stdx=1e-4
    stdy=1e-4

    kernelsize=3

    # get absolute peak pixel
    max_f = np.amax(f)
    [xpeak,ypeak] = np.unravel_index(f.argmax(), f.shape)

    if subpixel==False or xpeak < kernelsize or xpeak > np.shape(f)[0]-kernelsize or ypeak < kernelsize or ypeak > np.shape(f)[1]-kernelsize: # on edge
        return xpeak, ypeak, stdx, stdy, max_f # return absolute peak
    else:
        # determine sub pixel accuracy by centroid
        fextracted=f[xpeak-kernelsize:xpeak+kernelsize+1,ypeak-kernelsize:ypeak+kernelsize+1]
        fextractedx=np.mean(fextracted,0)
        fextractedy=np.mean(fextracted,1)
        x=np.arange(-kernelsize,kernelsize+1,1)
        y=np.transpose(x)

        xoffset=np.dot(x,fextractedx)
        yoffset=np.dot(y,fextractedy)

        # return only one-thousandths of a pixel precision
        xoffset = np.round(1000*xoffset)/1000
        yoffset = np.round(1000*yoffset)/1000
        xpeak=xpeak+xoffset
        ypeak=ypeak+yoffset

        # 2D linear interpolation
        bilinterp = interpolate.interp2d(x, x, fextracted, kind='linear')
        max_f = bilinterp(xoffset,yoffset)

        # peak width (full width at half maximum)
        scalehalfwidth=1.1774
        stdx=scalehalfwidth*np.std(fextractedx)
        stdy=scalehalfwidth*np.std(fextractedy)

    return xpeak, ypeak, stdx, stdy, max_f 

def surveillance_error(targets, predictions):
    data = np.loadtxt('seg.csv')

    xs = np.linspace(0, 600, data.shape[0])
    ys = np.linspace(0, 600, data.shape[1])
    f = interp2d(xs, ys, np.transpose(data))

    scores = np.concatenate([f(t, p) for (t, p) in zip(targets, predictions)])
    return np.mean(scores), np.std(scores) 

def dense_image_warp(image, flow):
    # batch_size, height, width, channels = (array_ops.shape(image)[0],
    #                                        array_ops.shape(image)[1],
    #                                        array_ops.shape(image)[2],
    #                                        array_ops.shape(image)[3])
    batch_size, height, width, channels = (np.shape(image)[0],
                                           np.shape(image)[1],
                                           np.shape(image)[2],
                                           np.shape(image)[3])

    # The flow is defined on the image grid. Turn the flow into a list of query
    # points in the grid space.
    # grid_x, grid_y = array_ops.meshgrid(
    #     math_ops.range(width), math_ops.range(height))
    # stacked_grid = math_ops.cast(
    #     array_ops.stack([grid_y, grid_x], axis=2), flow.dtype)
    # batched_grid = array_ops.expand_dims(stacked_grid, axis=0)
    # query_points_on_grid = batched_grid - flow
    # query_points_flattened = array_ops.reshape(query_points_on_grid,
    #                                            [batch_size, height * width, 2])
    grid_x, grid_y = np.meshgrid(
        np.range(width), np.range(height))
    stacked_grid = np.cast(
        np.stack([grid_y, grid_x], axis=2), flow.dtype)
    batched_grid = np.expand_dims(stacked_grid, axis=0)
    query_points_on_grid = batched_grid - flow
    query_points_flattened = np.reshape(query_points_on_grid,
                                        [batch_size, height * width, 2])
    # Compute values at the query points, then reshape the result back to the
    # image grid.
    interpolated = interp2d(image, query_points_flattened)
    interpolated = np.reshape(interpolated,
                              [batch_size, height, width, channels])
    return interpolated 

def poly2d(self):
        #check for monotonicity ?
        #fix this, interp needs increasing
        poly2d = interp2d(self.size, self.alpha, self.crit_table)
        return poly2d 

def eps_func(self):
        '''
        function: a function that when passed a `x` and `y` values,
            returns the permittivity profile of the structure,
            interpolating if necessary.
        '''
        interp_real = interpolate.interp2d(self.x, self.y, self.eps.real)
        interp_imag = interpolate.interp2d(self.x, self.y, self.eps.imag)
        interp = lambda x, y: interp_real(x, y) + 1.j*interp_imag(x, y)
        return interp 

def n_func(self):
        '''
        function: a function that when passed a `x` and `y` values,
            returns the refractive index profile of the structure,
            interpolating if necessary.
        '''
        return interpolate.interp2d(self.x, self.y, self.n) 

def get_translate(workdir=None):

    filename = os.path.join(workdir, "Vicalloy/Fe-Co-V_140922a_META_DATA.csv")
    compdata_f = pd.read_csv(filename, sep="\t").dropna()
    print compdata_f.head()
    x = compdata_f["Xnom (mm)"].values
    y = compdata_f["Ynom (mm)"].values
    Co_concentration = compdata_f["Co (at%)"].values
    Fe_concentration = compdata_f["Fe (at%)"].values
    V_concentration = compdata_f["V (at%)"].values
    method = "linear"
    # method = 'nearest'

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        Co_concI = interp2d(x, y, Co_concentration, kind=method)
        Fe_concI = interp2d(x, y, Fe_concentration, kind=method)
        V_concI = interp2d(x, y, V_concentration, kind=method)

    def translate(key):
        manip_z, manip_y = key
        sample_y = manip_z - 69.5
        sample_x = (manip_y + 8) * 2
        Co = Co_concI(sample_x, sample_y)[0] / 100.0
        Fe = Fe_concI(sample_x, sample_y)[0] / 100.0
        V = V_concI(sample_x, sample_y)[0] / 100.0
        return ("Fe{:.2f}Co{:.2f}V{:.2f}".format(Fe, Co, V), sample_x, sample_y)

    return translate 

def interpolate_zvals_at_xy(xy, topography, method='interp2d'):
        """
        Interpolates DEM values on a defined section

        Args:
            xy: x (EW) and y (NS) coordinates of the profile
            topography (:class:`gempy.core.grid_modules.topography.Topography`)
            method: interpolation method, 'interp2d' for cubic scipy.interpolate.interp2d
                                             'spline' for scipy.interpolate.RectBivariateSpline

        Returns:
            numpy.ndarray: z values, i.e. topography along the profile

        """

        xj = topography.values_2d[:, 0, 0]
        yj = topography.values_2d[0, :, 1]
        zj = topography.values_2d[:, :, 2]

        if method == 'interp2d':
            f = interpolate.interp2d(xj, yj, zj.T, kind='cubic')
            zi = f(xy[:, 0], xy[:, 1])
            if xy[:, 0][0] <= xy[:, 0][-1] and xy[:, 1][0] <= xy[:, 1][-1]:
                return np.diag(zi)
            else:
                return np.flipud(zi).diagonal()
        else:
            assert xy[:, 0][0] <= xy[:, 0][-1], 'The xy values of the first point must be smaller than second.' \
                                               'Please use interp2d as method argument. Will be fixed.'
            assert xy[:, 1][0] <= xy[:, 1][-1], 'The xy values of the first point must be smaller than second.' \
                                               'Please use interp2d as method argument. Will be fixed.'
            f = interpolate.RectBivariateSpline(xj, yj, zj)
            zi = f(xy[:, 0], xy[:, 1])
            return np.flipud(zi).diagonal() 

def interpolate2d(r,v):
    """Return a function that does 2d interpolation"""
    from scipy.interpolate import interp2d
    x,y = r[:,0],r[:,1] # data
    grid_x, grid_y = np.mgrid[np.min(x):np.max(x):100j,np.min(y):np.max(y):100j]
    from scipy.interpolate import griddata
    grid_z = griddata(r,v, (grid_x, grid_y), method='nearest')
    f = interp2d(grid_x, grid_y, grid_z, kind='linear')
    return lambda ri: f(ri[0],ri[1]) 

def linStackedAmplitude_alt1(data,profilePos,twtt,vVals,tVals,typefact):
    '''
    Calculates the linear stacked amplitudes for each two-way 
    travel time sample and the provided velocity range 
    by summing the pixels of the data that follow a line given 
    by the two-way travel time zero offset and the velocity.

    INPUT:
    data          data matrix whose columns contain the traces
    profilePos    along-profile coordinates of the traces
    twtt          two-way travel time values for the samples, in ns
    vVals         list of velocity values for which to calculate the
                  linear stacked amplitudes, in m/ns
    tVals         list of twtt zero-offsets for which to calculate
                  the linear stacked amplitudes, in ns
    typefact      factor for antenna separation depending if this is
                  for CMP (typefact=2) or WARR (typefact=1) data

    OUTPUT:
    linStAmp      matrix containing the linear stacked amplitudes
                  for the given data, tVals, and vVals
    '''
    linStAmp=np.zeros((len(tVals),len(vVals)))
    f = interp.interp2d(profilePos, twtt, data)        
    for vi in  tqdm(range(0,len(vVals))):
        for ti in range(0,len(tVals)):
            t = tVals[ti] + typefact*profilePos/vVals[vi]            
            vals = np.diagonal(np.asmatrix(f(profilePos, t)))
            linStAmp[ti,vi] = np.abs(sum(vals)/len(vals))
    return linStAmp 

def smoothNoise(noise, scale):
    """
    Given an 2D array of random values, it creates and array of
    same size. The values of the new arrays have lower variance
    between its neighbourghs. This is done by taking a subset of the array
    of size noise.size()/scale and extrapolating it to the original size
    
    Arguments:
        noise (array of float): A 2D array of floats to be smoothened
        scale (int): non negative. An extrapolation factor.
            e.g. scale = 2 will take quarter of the original image and 
            extrapolate it to the original size
        
    Returns:
        smoothed (array of float): A 2D array of same size as noise,
                                    with lower variance between negihbourghs
    
    """

    x = np.linspace(0,scale,noise.shape[1])
    y = np.linspace(0,scale,noise.shape[0])

    noise_scaled = noise[:len(x),:len(y)]

    X = np.linspace(0,1,noise.shape[1])
    Y = np.linspace(0,1,noise.shape[0])

    N = interp2d(x,y,noise_scaled)
    smoothed = N(X,Y)

    return smoothed 

def get_concentration_functions(composition_table_dict):

    meta = composition_table_dict["meta"]
    composition_table = Table.from_dict(composition_table_dict["data"])
    elements = [col for col in composition_table.columns if col not in meta]
    x = composition_table["X"].values
    y = composition_table["Y"].values
    cats = composition_table["X"].unique()
    concentration, conc, d, y_c, functions = {}, {}, {}, {}, RecursiveDict()

    for el in elements:
        concentration[el] = to_numeric(composition_table[el].values) / 100.0
        conc[el], d[el], y_c[el] = {}, {}, {}

        if meta["X"] == "category":
            for i in cats:
                k = "{:06.2f}".format(float(i))
                y_c[el][k] = to_numeric(y[where(x == i)])
                conc[el][k] = to_numeric(concentration[el][where(x == i)])
                d[el][k] = interp1d(y_c[el][k], conc[el][k])

            functions[el] = lambda a, b, el=el: d[el][a](b)

        else:
            functions[el] = interp2d(float(x), float(y), concentration[el])

    return functions 

def create_vtgt_const(self, v_tgt):
        self.vtgt[0].fill(v_tgt[0])
        self.vtgt[1].fill(v_tgt[1])        

        self.vtgt_interp_x = interpolate.interp2d(self.map[0,:,0], self.map[1,0,:], self.vtgt[0].T, kind='linear')
        self.vtgt_interp_y = interpolate.interp2d(self.map[0,:,0], self.map[1,0,:], self.vtgt[1].T, kind='linear') 

def __init__(self,orbit_datapath=None,f_nStars=10,**specs): 

        ObservatoryL2Halo.__init__(self,**specs)  
        self.f_nStars = int(f_nStars)
        
        # instantiating fake star catalog, used to generate good dVmap
        fTL = TargetList(**{"ntargs":self.f_nStars,'modules':{"StarCatalog": "FakeCatalog", \
                    "TargetList":" ","OpticalSystem": "Nemati", "ZodiacalLight": "Stark", "PostProcessing": " ", \
                    "Completeness": " ","BackgroundSources": "GalaxiesFaintStars", "PlanetPhysicalModel": " ", \
                    "PlanetPopulation": "KeplerLike1"}, "scienceInstruments": [{ "name": "imager"}],  \
                    "starlightSuppressionSystems": [{ "name": "HLC-565"}]   })
        
        f_sInds = np.arange(0,fTL.nStars)
        dV,ang,dt = self.generate_dVMap(fTL,0,f_sInds,self.equinox[0])
        
        # pick out unique angle values
        ang, unq = np.unique(ang, return_index=True)
        dV = dV[:,unq]
        
        #create dV 2D interpolant
        self.dV_interp  = interp.interp2d(dt,ang,dV.T,kind='linear') 

def create_vtgt_sink(self, p_sink, d_sink, v_amp_rng, v_phase0=np.random.uniform(-np.pi, np.pi)):
        # set vtgt orientations
        rng_xy = (-p_sink + self.map_rng_xy.T).T
        self.vtgt = -self._generate_grid(rng_xy, self.res_map)

        # set vtgt amplitudes
        self._set_sink_vtgt_amp(p_sink, d_sink, v_amp_rng, v_phase0)

        self.vtgt_interp_x = interpolate.interp2d(self.map[0,:,0], self.map[1,0,:], self.vtgt[0].T, kind='linear')
        self.vtgt_interp_y = interpolate.interp2d(self.map[0,:,0], self.map[1,0,:], self.vtgt[1].T, kind='linear')

# ----------------------------------------------------------------------------------------------------------------- 

def test_interp2d_linear(self):
        # Ticket #898
        a = np.zeros([5, 5])
        a[2, 2] = 1.0
        x = y = np.arange(5)
        b = interp2d(x, y, a, 'linear')
        assert_almost_equal(b(2.0, 1.5), np.array([0.5]), decimal=2)
        assert_almost_equal(b(2.0, 2.5), np.array([0.5]), decimal=2) 

def test_interp2d_meshgrid_input(self):
        # Ticket #703
        x = linspace(0, 2, 16)
        y = linspace(0, pi, 21)
        z = sin(x[None,:] + y[:,None]/2.)
        I = interp2d(x, y, z)
        assert_almost_equal(I(1.0, 2.0), sin(2.0), decimal=2) 

def plot_mcontour(self, ndim0, ndim1, z, show_mode):
        "use mayavi.mlab to plot contour."
        if not mayavi_installed:
            self.__logger.info("Mayavi is not installed on your device.")
            return
        #do 2d interpolation
        #get slice object
        s = np.s_[0:ndim0:1, 0:ndim1:1]
        x, y = np.ogrid[s]
        mx, my = np.mgrid[s]
        #use cubic 2d interpolation
        interpfunc = interp2d(x, y, z, kind='cubic')
        newx = np.linspace(0, ndim0, 600)
        newy = np.linspace(0, ndim1, 600)
        newz = interpfunc(newx, newy)
        #mlab
        face = mlab.surf(newx, newy, newz, warp_scale=2)
        mlab.axes(xlabel='x', ylabel='y', zlabel='z')
        mlab.outline(face)
        #save or show
        if show_mode == 'show':
            mlab.show()
        elif show_mode == 'save':
            mlab.savefig('mlab_contour3d.png')
        else:
            raise ValueError('Unrecognized show mode parameter : ' +
                             show_mode)

        return