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

def __init__(self, A, h, alpha, CL, CD, rho, smoothing=0, k_spline=3):
        self.A   = A
        self.h   = h
        self.Asp = 2*self.h**2/self.A

        self.rho = rho

        # Input lift and drag data
        self.n     = len(alpha)
        self.alpha = alpha
        self.CL    = CL
        self.CD    = CD

        self.k_spline = k_spline

        # -------- Spline interpolation -----------------------------
        if len(self.alpha.shape) == 1:
            self.interpolationMethod = 'spline'
            self.nrControlParameters = 1
            self.CL_spline = interpolate.splrep(self.alpha, self.CL, s=smoothing, k=self.k_spline)
            self.CD_spline = interpolate.splrep(self.alpha, self.CD, s=smoothing, k=self.k_spline)
        elif len(self.alpha.shape) == 2:
            self.interpolationMethod = 'griddata'
            self.nrControlParameters = 2
            self.CL_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CL, smooth=smoothing)
            self.CD_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CD, smooth=smoothing)

def CLCD(self, alpha):
        if self.interpolationMethod == 'spline':
            CL = interpolate.splev(alpha, self.CL_spline)
            CD = interpolate.splev(alpha, self.CD_spline)
        elif self.interpolationMethod == 'Rbf':
            CL = self.CL_RbfModel(alpha[0], alpha[1])
            CD = self.CD_RbfModel(alpha[0], alpha[1])
        elif self.interpolationMethod == 'griddata':
            CL = interpolate.griddata(self.alpha, self.CL, alpha)
            CD = interpolate.griddata(self.alpha, self.CD, alpha)

        if self.nrControlParameters == 2:
            return float(CL), float(CD)
        else:
            return CL, CD

def populateMissingData(self,approach="Smooth",ilog=None):
        '''
        This function is used to interpolate missing data in the image.
        '''
        if approach == 'Smooth':
            # first run a median filter over the array, then smooth the result.
            #xmin,xmax,ymin,ymax,zmin,zmax = self.getRange()
            mask = np.array(self.flags,dtype=np.bool)

            z = self.getZImage().asMatrix2D()
            median = nd.median_filter(z,size=(15,15))

            mask = mask.flatten()
            z = z.flatten()
            median = median.flatten()

            z[ mask==False ] = median[ mask==False ]

            if ilog != None:
                ilog.log(pv.Image(median.reshape(self.width,self.height)),label="Median")
                ilog.log(pv.Image(z.reshape(self.width,self.height)),label="ZMedian")

            mask = mask.flatten()
            z = z.flatten()
            median = median.flatten()

            for i in range(5):
                tmp = z.copy()
                smooth = nd.gaussian_filter(z.reshape(self.width,self.height),2.0).flatten()
                z[ mask==False ] = smooth[ mask==False ]
                print "Iteration:",i,(z-tmp).max(),(z-tmp).min()
                ilog.log(pv.Image(z.reshape(self.width,self.height)),label="ZSmooth%02d"%i)
                ilog.log(pv.Image((z-tmp).reshape(self.width,self.height)),label="ZSmooth%02d"%i)


        if approach == 'RBF':
            mask = np.array(self.flags,dtype=np.bool)
            mask = mask.flatten()

            x = np.arange(self.width).reshape(self.width,1)
            x = x*np.ones((1,self.height))
            x = x.flatten()

            y = np.arange(self.height).reshape(1,self.height)
            y = y*np.ones((self.width,1))
            y = y.flatten()

            z = self.z.copy()
            z = z.flatten()

            print "Coords:"
            print len(mask)
            print len(x[mask])
            print len(y[mask])
            print len(z[mask])

            # this produces an error.  Probably has too much data
            it.Rbf(x[mask],y[mask],z[mask])
            pass

def compute_ck2004_response(self, path, verbose=False):
        """
        Computes Castelli & Kurucz (2004) intensities across the entire
        range of model atmospheres.

        @path: path to the directory containing ck2004 SEDs
        @verbose: switch to determine whether computing progress should
        be printed on screen

        Returns: n/a
        """

        models = glob.glob(path+'/*M1.000*')
        Nmodels = len(models)

        # Store the length of the filename extensions for parsing:
        offset = len(models[0])-models[0].rfind('.')

        Teff, logg, abun = np.empty(Nmodels), np.empty(Nmodels), np.empty(Nmodels)
        InormE, InormP = np.empty(Nmodels), np.empty(Nmodels)

        if verbose:
            print('Computing Castelli & Kurucz (2004) passband intensities for %s:%s. This will take a while.' % (self.pbset, self.pbname))

        for i, model in enumerate(models):
            #~ spc = np.loadtxt(model).T -- waaay slower
            spc = np.fromfile(model, sep=' ').reshape(-1,2).T

            Teff[i] = float(model[-17-offset:-12-offset])
            logg[i] = float(model[-11-offset:-9-offset])/10
            sign = 1. if model[-9-offset]=='P' else -1.
            abun[i] = sign*float(model[-8-offset:-6-offset])/10

            spc[0] /= 1e10 # AA -> m
            spc[1] *= 1e7  # erg/s/cm^2/A -> W/m^3
            wl = spc[0][(spc[0] >= self.ptf_table['wl'][0]) & (spc[0] <= self.ptf_table['wl'][-1])]
            fl = spc[1][(spc[0] >= self.ptf_table['wl'][0]) & (spc[0] <= self.ptf_table['wl'][-1])]
            fl *= self.ptf(wl)
            flP = fl*wl
            InormE[i] = np.log10(fl.sum()/self.ptf_area*(wl[1]-wl[0]))             # energy-weighted intensity
            InormP[i] = np.log10(flP.sum()/self.ptf_photon_area*(wl[1]-wl[0]))     # photon-weighted intensity
            if verbose:
                if 100*i % (len(models)) == 0:
                    print('%d%% done.' % (100*i/(len(models)-1)))

        # Store axes (Teff, logg, abun) and the full grid of Inorm, with
        # nans where the grid isn't complete.
        self._ck2004_axes = (np.unique(Teff), np.unique(logg), np.unique(abun))

        self._ck2004_energy_grid = np.nan*np.ones((len(self._ck2004_axes[0]), len(self._ck2004_axes[1]), len(self._ck2004_axes[2]), 1))
        self._ck2004_photon_grid = np.nan*np.ones((len(self._ck2004_axes[0]), len(self._ck2004_axes[1]), len(self._ck2004_axes[2]), 1))
        for i, I0 in enumerate(InormE):
            self._ck2004_energy_grid[Teff[i] == self._ck2004_axes[0], logg[i] == self._ck2004_axes[1], abun[i] == self._ck2004_axes[2], 0] = I0
        for i, I0 in enumerate(InormP):
            self._ck2004_photon_grid[Teff[i] == self._ck2004_axes[0], logg[i] == self._ck2004_axes[1], abun[i] == self._ck2004_axes[2], 0] = I0

        # Tried radial basis functions but they were just terrible.
        #~ self._log10_Inorm_ck2004 = interpolate.Rbf(self._ck2004_Teff, self._ck2004_logg, self._ck2004_met, self._ck2004_Inorm, function='linear')
        self.content.append('ck2004')
        self.atmlist.append('ck2004')