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

def __call__(self, vals, fill_value=np.nan):
        """
        Evaluate interpolator for values given at the source points.

        Parameters
        ----------
        vals : ndarray of float, shape (numsourcepoints, ...)
            Values at the source points which to interpolate
        fill_value : float
            is needed if linear interpolation fails; defaults to np.nan

        Returns
        -------
        output : ndarray of float with shape (numtargetpoints,...)

        """
        self._check_shape(vals)
        ip = LinearNDInterpolator(self.src, vals, fill_value=fill_value)
        return ip(self.trg)


# -----------------------------------------------------------------------------
# Covariance routines needed for Kriging
# -----------------------------------------------------------------------------

def disparity_interpolate(disp, fill_zeros=True): 
    # Determine valid positive disparity pixels
    xyd = valid_pixels(disp, disp > 0)

    # Nearest neighbor interpolator
    nn = LinearNDInterpolator(xyd[:,:2], xyd[:,2])

    # Interpolate all pixels
    xyd = valid_pixels(disp, np.ones(shape=disp.shape, dtype=np.bool))
    return nn(xyd[:,:2]).reshape(disp.shape[:2])

def __init__(
            self,
            name,
            grid_size,
            x_label,
            y_label,
            x_lim,
            y_lim,
            file_name,
            order_parameter
    ):
        self.name = name
        self.grid_size = grid_size
        self.x_lim = x_lim
        self.y_lim = y_lim
        self.x_label = x_label
        self.y_label = y_label
        self.file_name = file_name
        self.order_parameter = order_parameter
        # load data file
        data = np.load(file_name)
        shape = (grid_size, grid_size)
        # for background plotting
        xx = np.reshape(data[:, 0], shape)
        yy = np.reshape(data[:, 1], shape)
        # color plot
        zz = [
            order_parameter(entry)
            for entry in data[:, 2:]
        ]
        zz = np.reshape(zz, shape)
        self.plot_data = xx, yy, zz
        # interpolation
        # rescale parameter space to (0,1)
        data[:, 0] -= np.min(data[:, 0])
        data[:, 0] /= np.max(data[:, 0])
        data[:, 1] -= np.min(data[:, 1])
        data[:, 1] /= np.max(data[:, 1])
        self._interp = LinearNDInterpolator(data[:, 0:2], data[:, 2:], fill_value=np.nan)

def create_interpolator(filename):
  """
  Loads a LUT file and creates an interpolated LUT object.
  The interpolant is constructed by triangulating the input data using Qhull
  and performing linear barycentric interpolation on each triangle

  """

  #load LUT
  LUT = pickle.load(open(filename,"rb"))

  # LUT inputs (H2O, O3, etc.) and outputs (i.e. atmcorr coeffs)
  inputs = permutate_invars(LUT['config']['invars'])
  outputs = LUT['outputs']

  # piecewise linear interpolant in N dimensions
  t = time.time()
  interpolator = LinearNDInterpolator(inputs,outputs)
  print('Interpolation took {:.2f} (secs) = '.format(time.time()-t))

  # sanity check
  print('Quick check..')
  i = 0
  true   = (outputs[i][0],outputs[i][1])
  interp = interpolator(inputs[i][0],inputs[i][1],inputs[i][2],inputs[i][3],inputs[i][4]) 
  print('true   = {0[0]:.2f} {0[1]:.2f}'.format(true))
  print('interp = {0[0]:.2f} {0[1]:.2f}'.format(interp))

  return interpolator

def bellman_operator(self, v):
        """

        The Bellman operator.  Including for comparison. Value function
        iteration is not recommended for this problem.  See the
        reservation wage operator below.

        Parameters
        ----------
        v : array_like(float, ndim=1, length=len(?_grid))
            An approximate value function represented as a
            one-dimensional array.

        Returns
        -------
        new_v : array_like(float, ndim=1, length=len(?_grid))
            The updated value function

        """
        # == Simplify names == #
        f, g, ?, c, q = self.f, self.g, self.?, self.c, self.q

        vf = LinearNDInterpolator(self.grid_points, v)
        N = len(v)
        new_v = np.empty(N)

        for i in range(N):
            w, ? = self.grid_points[i, :]
            v1 = w / (1 - ?)
            integrand = lambda m: vf(m, q(m, ?)) * (? * f(m)
                                                     + (1 - ?) * g(m))
            integral, error = fixed_quad(integrand, 0, self.w_max)
            v2 = c + ? * integral
            new_v[i] = max(v1, v2)

        return new_v

def get_greedy(self, v):
        """
        Compute optimal actions taking v as the value function.

        Parameters
        ----------
        v : array_like(float, ndim=1, length=len(?_grid))
            An approximate value function represented as a
            one-dimensional array.

        Returns
        -------
        policy : array_like(float, ndim=1, length=len(?_grid))
            The decision to accept or reject an offer where 1 indicates
            accept and 0 indicates reject

        """
        # == Simplify names == #
        f, g, ?, c, q = self.f, self.g, self.?, self.c, self.q

        vf = LinearNDInterpolator(self.grid_points, v)
        N = len(v)
        policy = np.zeros(N, dtype=int)

        for i in range(N):
            w, ? = self.grid_points[i, :]
            v1 = w / (1 - ?)
            integrand = lambda m: vf(m, q(m, ?)) * (? * f(m) +
                                                     (1 - ?) * g(m))
            integral, error = fixed_quad(integrand, 0, self.w_max)
            v2 = c + ? * integral
            policy[i] = v1 > v2  # Evaluates to 1 or 0

        return policy

def get_norm_for_gauss_power(siglow=1.,sighigh=4.,powlow=1.5,powhigh=4.):
    S,P = np.meshgrid(np.linspace(siglow,sighigh,31), 
                      np.linspace(powlow,powhigh,26))
    x = np.linspace(-20,20,41)

    S = S.ravel()
    P = P.ravel()
    v = np.zeros(S.shape)
    for i,s in enumerate(S):
        v[i] = gauss_power(x, s, P[i]).sum()
    X = np.zeros((len(S),2))
    X[:,0] = S
    X[:,1] = P
    return LinearNDInterpolator(X,v, fill_value=1.)

def interpolate_LUTs(self):
    """
    interpolate look up tables
    """

    filepaths = sorted(glob.glob(self.LUTs_dir+os.path.sep+'*.lut'))

    if filepaths:

      print('running n-dimensional interpolation may take a few minutes...')

      try:

        for fpath in filepaths:

          fname = os.path.basename(fpath)
          fid, ext = os.path.splitext(fname)
          ilut_filepath = os.path.join(self.iLUTs_dir,fid+'.ilut')

          if os.path.isfile(ilut_filepath):
            print('iLUT file already exists (skipping interpolation): {}'.format(fname))
          else:
            print('Interpolating: '+fname)

            # load look up table
            LUT = pickle.load(open(fpath,"rb"))

            # input variables (all permutations)
            invars = LUT['config']['invars']
            inputs = list(product(invars['solar_zs'],
                                  invars['H2Os'],
                                  invars['O3s'],
                                  invars['AOTs'],
                                  invars['alts']))  

            # output variables (6S correction coefficients)
            outputs = LUT['outputs']

            # piecewise linear interpolant in n-dimensions
            t = time.time()
            interpolator = LinearNDInterpolator(inputs,outputs)
            print('Interpolation took {:.2f} (secs) = '.format(time.time()-t))

            # save new interpolated LUT file
            pickle.dump(interpolator, open(ilut_filepath, 'wb' ))

      except:

        print('interpolation error')

    else:

      print('LUTs directory: ',self.LUTs_dir)
      print('LUT files (.lut) not found in LUTs directory, try downloading?')

def interpolate_LUTs(self):
    """
    interpolates look up table files (.lut)
    """

    filepaths = sorted(glob.glob(self.LUT_path+os.path.sep+'*.lut'))
    if filepaths:
      print('\n...Running n-dimensional interpolation may take a several minutes (only need to do this once)...')
      try:
        for fpath in filepaths:
          fname = os.path.basename(fpath)
          fid, ext = os.path.splitext(fname)
          ilut_filepath = os.path.join(self.iLUT_path,fid+'.ilut')
          if os.path.isfile(ilut_filepath):
            print('iLUT file already exists (skipping interpolation): {}'.format(fname))
          else:
            print('Interpolating: '+fname)

            # load look up table
            LUT = pickle.load(open(fpath,"rb"))

            # input variables (all permutations)
            invars = LUT['config']['invars']
            inputs = list(product(invars['solar_zs'],
                                  invars['H2Os'],
                                  invars['O3s'],
                                  invars['AOTs'],
                                  invars['alts']))  

            # output variables (6S correction coefficients)
            outputs = LUT['outputs']

            # piecewise linear interpolant in n-dimensions
            t = time.time()
            interpolator = LinearNDInterpolator(inputs,outputs)
            print('Interpolation took {:.2f} (secs) = '.format(time.time()-t))

            # save new interpolated LUT file
            pickle.dump(interpolator, open(ilut_filepath, 'wb' ))

      except:

        print('interpolation error')

    else:

      print('look up tables files (.lut) not found in LUTs directory:\n{}'.format(self.LUT_path))