Python scipy.interpolate 模块，interp2d() 实例源码

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 _interp2_r0(Data, Pow2factor, kind, disp=False):
    if disp:
        p30 = np.poly1d(np.polyfit([10, 30, 50, 60, 70, 80],
                                   np.log([0.8, 1.3, 6, 12, 28, 56]), 2))
        print('Expectd time for NxN data:', np.exp(p30(Data.shape[0])))

    x = np.arange(Data.shape[1])
    y = np.arange(Data.shape[0])
    xv, yv = np.meshgrid(x, y)
    f = interpolate.interp2d(xv, yv, Data, kind=kind)

    xnew = np.arange(0, Data.shape[1], 1 / (2**Pow2factor))
    ynew = np.arange(0, Data.shape[0], 1 / (2**Pow2factor))
    Upsampled = f(xnew, ynew)

    return Upsampled

def eps_func(self):
        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):
        return interpolate.interp2d(self.x, self.y, self.n)

def __init__(self, today, proj, disc, factory, vols, mode):
        self.proj = proj
        self.disc = disc
        self.factory = factory
        self.today = today
        self.spotdate = factory.roll_date(today, 'spot') # also for marking startdate of hw vol term structure         
        self.surface = interp2d(vols.columns, vols.index, vols.values, kind='linear')
        self.mode = 'capfloor'    

        class Caplet(BlackLognormalModel if mode == 'lognormal' else BlackNormalModel):
            def __init__(self, libor, forward, blackvol, opt_term, annuity):
                self.underlier = libor 
                super(Caplet, self).__init__(forward, blackvol, opt_term, annuity)
        self.Caplet = Caplet

def bilin_interp(grid2d, ith, ir, dth, dr, nth, nr):
   '''
   This function finds an interpolated value in the 2d plume grid, and returns
   the value. It's much faster than scipy.interpolate.interp2d

   Takes as arguments the 2d grid of plume perturbations, the index in theta,
   the index in radius, step size in theta and radius, and number of points
   in theta and radius.
   '''
   value = grid2d[ith,ir]*(1-dth)*(1-dr) + grid2d[ith+1,ir]*dth*(1-dr)
   value += grid2d[ith,ir+1]*(1-dth)*dr + grid2d[ith+1,ir+1]*dth*dr
   return value

def __init__(self, input_data, output_data, name=""):
        # check type and dimensions
        assert isinstance(input_data, list)
        assert isinstance(output_data, np.ndarray)

        # output_data has to contain at least len(input_data) dimensions
        assert len(input_data) <= len(output_data.shape)

        for dim in range(len(input_data)):
            assert len(input_data[dim]) == output_data.shape[dim]

        self.input_data = input_data
        if output_data.size == 0:
            raise ValueError("No initialisation possible with an empty array!")
        self.output_data = output_data
        self.min = output_data.min()
        self.max = output_data.max()
        self.name = name

        # self._interpolator = si.interp2d(input_data[0], input_data[1], output_data, bounds_error=True)

def __init__(self, eval_data_list, time_point=None, spatial_point=None, ylabel="",
                 legend_label=None, legend_location=1, figure_size=(10, 6)):

        if not ((isinstance(time_point, Number) ^ isinstance(spatial_point, Number)) and \
                    (isinstance(time_point, type(None)) ^ isinstance(spatial_point, type(None)))):
            raise TypeError("Only one kwarg *_point can be passed,"
                            "which has to be an instance from type numbers.Number")

        DataPlot.__init__(self, eval_data_list)

        plt.figure(facecolor='white', figsize=figure_size)
        plt.ylabel(ylabel)
        plt.grid(True)

        # TODO: move to ut.EvalData
        len_data = len(self._data)
        interp_funcs = [si.interp2d(eval_data.input_data[1], eval_data.input_data[0], eval_data.output_data)
                        for eval_data in eval_data_list]

        if time_point is None:
            slice_input = [data_set.input_data[0] for data_set in self._data]
            slice_data = [interp_funcs[i](spatial_point, slice_input[i]) for i in range(len_data)]
            plt.xlabel('$t$')
        elif spatial_point is None:
            slice_input = [data_set.input_data[1] for data_set in self._data]
            slice_data = [interp_funcs[i](slice_input[i], time_point) for i in range(len_data)]
            plt.xlabel('$z$')
        else:
            raise TypeError

        if legend_label is None:
            show_leg = False
            legend_label = [evald.name for evald in eval_data_list]
        else:
            show_leg = True

        for i in range(0, len_data):
            plt.plot(slice_input[i], slice_data[i], label=legend_label[i])

        if show_leg:
            plt.legend(loc=legend_location)

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 barycentric_trapezoidial_interpolation(Nx,Ny,p,hexagonalOffset=0.5):
    '''
    define our function to calculate the position of points from Nx columns and Ny rows.
    The vertices are defined by p which is a size(4,2) array.
    each row of p are the coordinates or the vertices of our trapezoid
    the vertices have to be given in a specific order:
    [[1 1]
     [1 2]
     [2 1]
     [2 2]]
    an example plot using the barycentric interpolation to regrid data
    define number of rows and number of columns and the vertices, then make some plots

    Example:
    Nx = 20
    Ny = 15

    coords = [[0,0],[0,1],[1,0],[1,1]] #these are the [x,y] coords of your 4 draggable corners
    coords = np.asarray(coords)

    f, ax = plt.subplots(2, 2) # sharey=True, sharex=True)
    for i,a in enumerate(ax.flatten()):
        newCoords = coords[:]
        if i > 0:
            newCoords = newCoords + np.random.rand(4,2) / 5
        xi,yi = barycentric_trapezoidial_interpolation(Nx,Ny,newCoords)
        a.plot(xi,yi,'.',markersize=12)
    plt.show()
    '''     
    x_basis = np.linspace(0,1,Nx)
    y_basis = np.linspace(0,1,Ny)

    px = [[p[0,0], p[2,0]],[p[1,0], p[3,0]]] #these are the [2,2] x-coordinates
    py = [[p[0,1], p[2,1]],[p[1,1], p[3,1]]] #these are the [2,2] y-coordinates
    fx = interpolate.interp2d([0,1], [0,1], px, kind='linear')
    xi = fx(x_basis[:],y_basis[:]).flatten()
    fy = interpolate.interp2d([0,1], [0,1], py, kind='linear')
    yi = fy(x_basis[:],y_basis[:]).flatten()
    d1 = (p[2,0] - p[0,0]) / Nx / 2.0
    d2 = (p[3,0] - p[1,0]) / Nx / 2.0
    offset = (d1 + d2) * hexagonalOffset
    #every other row will be shifted in diff(x) * hexagonalOffset
    for i in range(0,len(xi)-Nx,Nx*2):
        for j in range(Nx):
            xi[i+j+Nx] += offset
    return xi,yi

def interpolate_curvature(curvatures, times, kind='cubic'):
    """ Interpolate the curvatures.

    A curvature is a parametrisation `s`, d(angle)/ ds and a parameter between [0,1].
    Return a surface f(s,t).
    """
    import numpy as np

    from scipy import interpolate
    from scipy.ndimage import measurements

    if len(curvatures) == 3 and not isinstance(curvatures[2], tuple):
        curvatures= [curvatures]

    curv_abs = [c[0] for c in curvatures]
    curves = [c[1] for c in curvatures]
    params = [c[2] for c in curvatures]

    n = len(params)
    if n==1:
        curv_abs*=2
        curves*=2
        curves[0]=np.zeros(len(curves[0]))
        params = [0.,1.]
    elif False:
        if params[0] != 0.:
            params.insert(0,0.)
            curves.insert(0,curves[0])
            curv_abs.insert(0, curv_abs[0])
        if params[-1] != 0:
            params.append(1.)
            curves.append(curves[-1])
            curv_abs.append(curv_abs[-1])

    # compute a common parametrisation s
    # We conserve the last parametrisation because the result is very sensitive
    if True:
        min_s = min(np.diff(s).min() for s in curv_abs)
        s_new = np.unique(np.array(curv_abs).flatten())
        ds = np.diff(s_new)
        k = np.cumsum(ds >= min_s)
        labels = range(k.min(), k.max()+1)
        s = np.zeros(len(labels)+1)
        s[1:] = measurements.mean(s_new[1:], k, labels)
        s[-1]=1.
        s = s_new
    else:
        s = np.array(curv_abs[-1])

    # renormalise all the curves
    curves = [np.interp(s[1:-1], old_s[1:-1], old_crv) for old_s, old_crv in zip(curv_abs,curves)]

    # interpolate correctly the curvatures
    x = s[1:-1]
    y = np.array(params)
    z = np.array(curves)

    #f = interpolate.interp2d(y, x, z, kind=kind)
    f = interpolate.RectBivariateSpline(x, y, z.T, kx=1, ky=1)
    return s, f(x,times)