Python interpolate.bisplev函数代码示例

 def __interpolateParameters__(self, height, latitude, tkcDeep, tkcUpper):
     # Preallocate result array and start looping through all values
     isScalar = not util.isArray(height)
     results = []
     
     if isScalar:
         results = [0]
         height = [height]
         latitude = [latitude]
     else:
         results = np.zeros(height.shape)
     
     for i in range(0, len(height)):
         # Check where the height is with respect to the interpolation limits
         if height[i] <= tkcDeep[0][-1]:
             results[i] = scp_ip.splev(height[i], tkcDeep)
         elif height[i] >= tkcUpper[0][0]:
             results[i] = scp_ip.bisplev(height[i], latitude[i], tkcUpper)
         else:
             # Interpolate between the lower and upper interpolating functions (do
             # so linearly for now)
             low = scp_ip.splev(tkcDeep[0][-1], tkcDeep)
             high = scp_ip.bisplev(tkcUpper[0][0], latitude[i], tkcUpper)
             
             results[i] = low + (high - low) * (height[i] - tkcDeep[0][-1]) / \
                 (tkcUpper[0][0] - tkcDeep[0][-1])
                 
     if isScalar:
         return results[0]
         
     return results

    def apar_intrp(x,y,z):

        dx = dx_xy[x,y]
        dy = dy_xy[x,y]
        
        
        # Interpolating down the y-axis (along the field)
        y_vals = np.array(range(ny), dtype=float)

        for j in range(len(y_vals)):
            y_vals[j] = interpolate.bisplev(x,z,tck[j])
#        print y_vals, 'sdf'
        dy_coeffs = interpolate.splrep(range(ny), y_vals, k=3)
        
        # Interpolating along the slices of data
#        intrp = interpolate.RectBivariateSpline(range(nx),range(nz),data[:,y,:],kx=3,ky=3)
        
#        tx,ty = intrp.get_knots()
#        tck = (tx,ty,intrp.get_coeffs(),3,3)
#        print y, int(y), np.shape(data[:,y,:])
        # From the cubic spline coefficients, returns derivatives
#        print 's', int(np.rint(y)), int(y)
        dervs = ( interpolate.bisplev(x,z,tck[int(np.rint(y))], dx=1, dy=0)/dx,
                  interpolate.splev(y,dy_coeffs,der=1)/dy,
                  interpolate.bisplev(x,z,tck[int(np.rint(y))], dx=0, dy=1)/dz )

        return  dervs

    def make_pixel_lut(self, dims):
        """
        Generate an x and y image which maps the array indices into
        floating point array indices (to be corrected for pixel size later)

        returns 
        FIXME - check they are the right way around
                add some sort of known splinefile testcase
        """
        # Cache the value in case of multiple calls
        if self.pixel_lut is None:
            x_im = numpy.outer(range(dims[0]), numpy.ones(dims[1]))
            y_im = numpy.outer(numpy.ones(dims[0]), range(dims[1]))
            # xcor is tck2
            x_im = numpy.add( x_im,
                              bisplev( range(dims[1]),
                                               range(dims[0]),
                                               self.tck2 ).T,
                              x_im)
            # ycor is tck1
            y_im = numpy.add( y_im,
                              bisplev( range(dims[1]),
                                               range(dims[0]),
                                               self.tck1 ).T,
                              y_im)
            self.pixel_lut = x_im, y_im
        return self.pixel_lut

def transferFactorCalculator (rho, pt):

    # Scipy bi-linear spline representation data object.
    # Cf. [https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.bisplrep.html]
    tck = [ np.array([    1.25,   1.25,   1.25,   1.25,    6.75,    6.75,    6.75,    6.75]), # Knots, x
            np.array([  450.,   450.,   450.,   450.,   1950.,   1950.,   1950.,   1950.]),   # Knots, y
            np.array([  0.18261346,  0.2174425 ,  0.37762574,  0.10284952,  0.26108651, # Spline coefficients
                       -0.14541588,  0.05701827,  0.27361263,  0.54531852,  0.77814774,
                        0.2479843 ,  0.23509468, -0.04597834, -0.47218929,  0.01928886,
                       -0.09066243]), 
            3, # Spline degree, x
            3] # Spline degree, y

    # Limits for the transfer factor map. Return zero if outside.
    limits = { 'rho': (  1.,    7.),
               'pt':  (400., 2000.) }

            # Check limits.
    if (limits['rho'][0] <= rho <= limits['rho'][1]) and \
       (limits['pt'] [0] <= pt  <= limits['pt'] [1]):
        # Calculate and return transfer factor.
        return interpolate.bisplev(rho, pt, tck)

    # Return fallback value.
    return 0.

 def fit(self, data, poldegree, swidth, sheight, threshold):
     if int(threshold) == -1:
         threshold = (int(data.mean()) * 10) / 7
     dims = data.shape
     xList = []
     yList = []
     zList = []
     for y in xrange(0, dims[0] - 1, sheight):
         for x in xrange(0, dims[1] - 1, swidth):
             view = data[y:y + sheight, x:x + swidth]
             flatIndex = numpy.argmax(view)
             yIdx, xIdx = numpy.unravel_index(flatIndex, view.shape)
             zValue = view[yIdx, xIdx]
             if zValue <= threshold:
                 xList.append(x + xIdx)
                 yList.append(y + yIdx)
                 zList.append(zValue)
     if len(xList) < (poldegree + 1) * (poldegree + 1):
         raise ValueError("Not enough reference points.")
     tck = interpolate.bisplrep(yList, xList, zList,
                                kx=poldegree, ky=poldegree,
                                xb=0, yb=0,
                                xe=int(dims[0]), ye=int(dims[1]))
     clipmin, clipmax = data.min(), threshold
     return interpolate.bisplev(range(dims[0]), range(dims[1]),
                                tck).clip(clipmin, clipmax)

    def interpolate(self, xi, yi):
        from scipy import interpolate
        # Need to write a norm function that calculates distance from a rib...

        """
        def interp(x1,x2,x3, x1i, x2i):
            spline = interpolate.Rbf(x1, x2, x3, function='thin-plate', smooth=0)
            return spline(x1i,x2i)

        try:
            zi = interp(self.points.x, self.points.y, self.points.z, xi, yi)
        except np.linalg.linalg.LinAlgError:
            zi = interp(self.points.y, self.points.x, self.points.z, yi, xi)
        """

        # Segfaults... Problems with the way scipy is compiled?
        tck = interpolate.bisplrep(self.points.x, self.points.y, self.points.z)
        zi = interpolate.bisplev(yi, xi, tck)


        """
        spline = interpolate.Rbf(self.points.x, self.points.y, self.points.z,
                                 function='thin-plate', smooth=0)
        zi = spline(xi,yi)
        """

        return zi

 def queryDepth(self, xq, yq):
     #return 1 data point
     depth = bisplev(xq, yq, self.tck)
     if np.isnan(depth):
         print "NaN returned for depth"
     else:
         return max(self.minz,min(depth, self.maxz))

def test_scipy_approx():
    """
    Test SciPy approximation of B-Spline surface
    :return: None
    """
    terrain_data = [
        [0.0, 0.0, 0.0], [0.0, 0.5, 0.4], [0.0, 1.0, 0.0],
        [0.5, 0.0, 0.2], [0.5, 0.5, 0.8], [0.5, 1.0, 0.2],
        [1.0, 0.0, 0.0], [1.0, 0.5, 0.4], [1.0, 1.0, 0.0]]
    tX = [item[0] for item in terrain_data]
    tY = [item[1] for item in terrain_data]
    tZ = [item[2] for item in terrain_data]
    from scipy import interpolate
    print('SciPy approximation ...')
    start_time = time.time()
    tck, fp, ior, msg = interpolate.bisplrep(tX, tY, tZ, kx=2, ky=2, full_output=1)
    end_time = time.time()
    print('Computed in {0} seconds.'.format(end_time - start_time))
    occ_bspline = convert.bspline.scipy_to_occ(tck)
    # Compute difference between original terrain data and B-Spline surface
    u_num = v_num = 50
    points = [[float(i)/u_num, float(j)/u_num, 0.0] for i in range(v_num+1) for j in range(u_num+1)]
    points = [(it[0], it[1], interpolate.bisplev(it[0], it[1], tck)) for it in points]
    # points = terrain_data

    display_results(occ_bspline, points)

 def evaluatePartialDerivativeV(self, x, y):
     result = np.empty(self.coeffElems)
     
     for i in range(self.coeffElems):
         result[i] = interpolate.bisplev(x, y, self.tcks[i], dy=1)
         
     return result

 def evaluate(self, x, y):
     result = np.empty(self.coeffElems)
     
     for i in range(self.coeffElems):
         result[i] = interpolate.bisplev(x, y, self.tcks[i])
         
     return result

    def derivatives(self, alpha, Re):

        # note: direct call to bisplev will be unnecessary with latest scipy update (add derivative method)
        tck_cl = self.cl_spline.tck[:3] + self.cl_spline.degrees  # concatenate lists
        tck_cd = self.cd_spline.tck[:3] + self.cd_spline.degrees

        dcl_dalpha = bisplev(alpha, Re, tck_cl, dx=1, dy=0)
        dcd_dalpha = bisplev(alpha, Re, tck_cd, dx=1, dy=0)

        if self.one_Re:
            dcl_dRe = 0.0
            dcd_dRe = 0.0
        else:
            dcl_dRe = bisplev(alpha, Re, tck_cl, dx=0, dy=1)
            dcd_dRe = bisplev(alpha, Re, tck_cd, dx=0, dy=1)

        return dcl_dalpha, dcl_dRe, dcd_dalpha, dcd_dRe

def loggtracks(masslimit,location,fileroot, metalstr, plot=True):



    #This is an array of all the masses that are part of the filenames
    massarrstr=['120.0', '85.0', '60.0', '40.0', '25.0', '20.0', '15.0', '12.0', '9.0', '7.0', '5.0', '4.0', '3.0', '2.5', '2.0', '1.7', '1.5', '1.25', '1.0', '0.9']



    #this is for each mass, read in the file and plot the evolutionary tracks
    for imass in range(masslimit):

        filemass = massarrstr[imass]

        filename = location+fileroot+filemass

        DataIn = np.genfromtxt(filename, dtype="float", unpack=True)
    
        time = DataIn[3,:]
        timesec = time*365*24*60*60
        Teff = 10**DataIn[ 6,:]
        logTeff = [math.log10(jj) for jj in Teff]
        L = 10**DataIn[ 5,:]*Lsun
        Rsquared = (L/(4*math.pi*sigma*Teff**4))
        Mass = DataIn[4,:]*Msun
        surfaceGrav = Mass*G/(Rsquared)
        logsurfaceGrav = [math.log10(ii) for ii in surfaceGrav]


        if plot == True:

            fadeplot(Teff, logsurfaceGrav, time)
        
        q0s = np.zeros(len(Teff))
        for timestep in range(len(Teff)):
            q0s[timestep] = interpolate.bisplev(Teff[timestep],logsurfaceGrav[timestep],gridq0s)
            
        
        
        for kk in range(len(q0s)):
            if q0s[kk] < 0:
                q0s[kk] = -q0s[kk]
                q0s[kk] = np.log10(q0s[kk])


        totalQ0s = np.trapz(q0s,timesec)
        
        logtotalQ0s = np.log10(totalQ0s)
        #print "Mass of: " + massarr[imass]+" Produces "+str(totalQ0s)+" Photons"
        logq0ints[imass] = totalQ0s
        #plt.plot(timesec/timesec[-1], q0s, 'k')
        titlestr= "Z= "+metalstr
        #plt.title(titlestr)
        #plt.xlabel("Stellar Lifetime")
        #plt.ylabel("Photons / Second")
        #plt.ylim([0,2e50])
    #plt.show()    
    return;

    def aproximate_terrain(self):
        """
        Try to aproximate terrain with bspline surface
        """

        tck,fp,ior,msg = interpolate.bisplrep(self.tX, self.tY, self.tZ, kx=5, ky=5, full_output=1)
        self.tck[(self.min_x, self.min_y, self.max_x, self.max_y)] = tck
        # Compute difference between original terrain data and b-spline surface
        self.tW = [abs(it[2] - interpolate.bisplev(it[0], it[1], tck)) for it in self.terrain_data]

def interpgrid(x,y, xlist,ylist, xmap, ymap, kx=3, ky=3, s=50):
   ''' for position x,y and a 2-D mapping map(list),
       i.e., xmap[xlist,ylist],ymap[xlist,ylist] given on a grid xlist,ylist; 
       the nearest xlist, ylist positions to each x,y pair are found and 
       interpolated to yield  mapx(x,y),mapy(x,y)
         
   x,y : rank-1 arrays of data points
   xlist, ylist, xmap, ymap: rank-1 arrays of data points
   
   +
   2008-08-24 NPMK (MSSL)
   '''
   from scipy import interpolate
   # check if the input is right data type
   # ... TBD
   
   # compute the Bivariate-spline coefficients
   # kx = ky =  3 # cubic splines (smoothing)
   task = 0 # find spline for given smoothing factor
   # s = 50 # spline goes through the given points
   # eps = 1.0e-6  (0 < eps < 1)
   
   #(tck_x, ems1) 
   tck_x = interpolate.bisplrep(xlist,ylist,xmap,kx=kx,ky=ky,s=s)
   #(fp1, ier1, msg1) = ems1
   #if ier1 in [1,2,3]: 
   #   print 'an error occurred computing the bivariate spline (xmap) '
   #   print ier1, msg1
   #   # raise error
   #   return None
   tck_y = interpolate.bisplrep(xlist,ylist,ymap,kx=kx,ky=ky,s=s) 
   #(fp2, ier2, msg2) = ems2
   #if ier2 in [1,2,3]: 
   #   print 'an error occurred computing the bivariate spline (ymap) '
   #   print ier2, msg2
   #   # raise error
   #   return None
   # compute the spline    
   
   xval = interpolate.bisplev(x, y, tck_x)
   yval = interpolate.bisplev(x, y, tck_y)
   
   return xval,yval

def sample(sgridxy, im, dd=(0,0), hscaling=1. ):
    xs = sgridxy[0].flatten()
    ys = sgridxy[1].flatten()
    assert(xs.size == ys.size)
    v = np.zeros(xs.size)
    for i in np.arange(v.size):
        v[i] = (1./(hscaling**np.sum(dd)))*flt(interpolate.bisplev(xs[i],ys[i], im, dd[0], dd[1]) )

    #return v.reshape(sgridxy[0].shape)
    return v