Python numpy.polyval函数代码示例

    def calibrate(self, title, rng=None, cal=None):
        """ Returns calibrated bin centers for saved histogram with given title.
        If the cal argument is given, title must exactly match a key in the file
        self.ofile.

        "title" title of histogram. If cal is not None then title must exactly
            match a key in the hdf5 file self.ofile.
        "rng" is the index range of elements to return
        "cal" is a set of polynomial coefficients for np.polyval. A linear
            scaling by a factor s would require cal=[s,0]
        """
        rng = rng or [None, None]
        if cal is None:
            bins = self.bcent(title, rng)
            return np.polyval(self.cal, bins)
        else:
            with h5py.File(self.ofile, 'r') as ofile:
                if title not in ofile.keys():
                    raise ValueError(title + ' is not in file ' + self.ofile)
                elif title[0] == 'b':
                    bins = self.bcent(title[1:], rng)
                    return np.polyval(cal, bins)
                else:
                    axis = ofile[title][rng[0]:rng[1]]
                    return np.polyval(cal, axis)

def showExamplePolyFit(xs,ys,fitDegree1 = 1,fitDegree2 = 2):
    pylab.figure()    
    pylab.plot(xs,ys,'r.',ms=2.0,label = "measured")

    # poly fit to noise
    coeeff = numpy.polyfit(xs, ys, fitDegree1)

    # Predict the curve
    pys = numpy.polyval(numpy.poly1d(coeeff), xs)

    se = mse(ys, pys)
    r2 = rSquared(ys, pys)

    pylab.plot(xs,pys, 'g--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree1,se,r2))

    # Poly fit to noise
    coeeffs = numpy.polyfit(xs, ys, fitDegree2)

    # Predict the curve
    pys = numpy.polyval(numpy.poly1d(coeeffs), xs)

    se = mse(ys, pys)
    r2 = rSquared(ys, pys)

    pylab.plot(xs,pys, 'b--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree2,se,r2))

    pylab.legend()

def cd_sphere_vector(Re):
    "Computes the drag coefficient of a sphere as a function of the Reynolds number Re."
    # Curve fitted after fig . A -56 in Evett & Liu :% " Fluid Mechanics & Hydraulics ",
    # Schaum ' s Solved Problems McGraw - Hill 1989.

    from numpy import log10,array,polyval
    CD = zeros_like(Re)
   
    CD = where(Re<0,0.0,0.0)     # condition 1
    
    CD = where((Re > 0.0) & (Re <=0.5),24/Re,CD) # condition 2

    p = array([4.22,-14.05,34.87,0.658])
    CD = where((Re > 0.5) & (Re <=100.0),polyval(p,1.0/Re),CD) #condition 3

    p = array([-30.41,43.72,-17.08,2.41])
    CD = where((Re >100.0)  & (Re <=1.0e4) ,polyval(p,1.0/log10(Re)),CD) #condition 4

    p = array([-0.1584,2.031,-8.472,11.932])
    CD = where((Re > 1.0e4)  &  (Re <=3.35e5),polyval(p,log10(Re)),CD) #condition 5

    CD = where((Re > 3.35e5) & (Re <=5.0e5),91.08*(log10(Re/4.5e5))**4 + 0.0764,CD) #condition 6

    p  = array([-0.06338,1.1905,-7.332,14.93])
    CD = where((Re > 5.05e5)  &  (Re <=8.0e6),polyval(p,log10(Re)),CD) #condition 7
    
    CD = where(Re>8.0e6,0.2,CD)  # condition 8

    return CD

    def _make_axes(self):
        '''Construct axes from calibration fields in header file
        '''
        xcalib = self.header.xcalibration
        ycalib = self.header.ycalibration

        xcalib_valid = struct.unpack('?', xcalib.calib_valid)

        if xcalib_valid:
            xcalib_order, = struct.unpack('>B', xcalib.polynom_order) # polynomial order
            px = xcalib.polynom_coeff[:xcalib_order+1]
            px = np.array(px[::-1]) # reverse coefficients to use numpy polyval
            pixels = np.arange(1, self.header.xdim + 1)
            px = np.polyval(px, pixels)
        else:
            px = np.arange(1, self.header.xdim + 1)

        ycalib_valid = struct.unpack('?', ycalib.calib_valid)

        if ycalib_valid:
            ycalib_order, = struct.unpack('>B', ycalib.polynom_order) # polynomial order
            py = ycalib.polynom_coeff[:ycalib_order+1]
            py = np.array(py[::-1]) # reverse coefficients to use numpy polyval
            pixels = np.arange(1, self.header.ydim + 1)
            py = np.polyval(py, pixels)
        else:
            py = np.arange(1, self.header.ydim + 1)

        self._xaxis = px
        self._yaxis = py

        return px, py

 def sampleR3(averagedist,boxdims):
     """low-discrepancy sampling using primes.
     The samples are evenly distributed with an average distance of averagedist inside the box with dimensions boxdims.
     Algorithim from "Geometric Discrepancy: An Illustrated Guide" by Jiri Matousek"""
     minaxis = numpy.argmin(boxdims)
     maxaxis = numpy.argmax(boxdims)
     meddimdist = numpy.sort(boxdims)[1]
     # convert average distance to number of samples.... do simple 3rd degree polynomial fitting...
     x = meddimdist/averagedist
     if x < 25.6:
         N = int(numpy.polyval([ -3.50181522e-01,   2.70202333e+01,  -3.10449514e+02, 1.07887093e+03],x))
     elif x < 36.8:
         N = int(numpy.polyval([  4.39770585e-03,   1.10961031e+01,  -1.40066591e+02, 1.24563464e+03],x))
     else:
         N = int(numpy.polyval([5.60147111e-01,  -8.77459988e+01,   7.34286834e+03, -1.67779452e+05],x))
     pts = numpy.zeros((N,3))
     pts[:,0] = numpy.linspace(0.0,meddimdist,N)
     pts[:,1] = meddimdist*numpy.mod(0.5+0.5*numpy.sqrt(numpy.arange(0,5.0*N,5.0)),1.0)
     pts[:,2] = meddimdist*numpy.mod(0.5+3*numpy.sqrt(numpy.arange(0,13.0*N,13.0)),1.0)
     if boxdims[minaxis] < meddimdist:
         pts = pts[pts[:,minaxis]<=boxdims[minaxis],:]
     if boxdims[maxaxis] > meddimdist:
         # have to copy across the max dimension
         numfullcopies = numpy.floor(boxdims[maxaxis]/meddimdist)
         oldpts = pts
         pts = numpy.array(oldpts)
         for i in range(int(numfullcopies)-1):
             oldpts[:,maxaxis] += meddimdist
             pts = numpy.r_[pts,oldpts]
         if boxdims[maxaxis]/meddimdist > numfullcopies:
             oldpts[:,maxaxis] += meddimdist
             pts = numpy.r_[pts,oldpts[oldpts[:,maxaxis]<=boxdims[maxaxis],:]]
     return pts

def corrNonlinGetPar(linearDet,nonLinearDet,order=2,data_0=0,
    correct_0=0,plot=False,returnCorrectedDet=False):
  """ Find parameters for non linear correction
    *linearDet* should be an 1D array of the detector that is linear
    *nonLinearDet* is the detector that is sussposed to be none linear
    *data_0" is an offset to use for the data (used only if plotting)"
    *correct_0* offset of the "linear detector"""
  p =  np.polyfit(nonLinearDet,linearDet,order)
  p[-1] = p[-1]-correct_0
  if plot:
    d = corrNonlin(nonLinearDet,p,data_0=data_0,correct_0=correct_0)
    plt.plot(linearDet,nonLinearDet,".",label="before correction")
    plt.plot(linearDet,d,".",label="after correction")
    poly_lin = np.polyfit(linearDet,d,1)
    xmin = min(linearDet.min(),0)
    xtemp = np.asarray( (xmin,linearDet.max()) )
    plt.plot(xtemp,np.polyval(poly_lin,xtemp),label="linear fit")
    plt.plot(linearDet,d-np.polyval(poly_lin,linearDet),
       ".",label="difference after-linear")
    plt.xlabel("linearDet")
    plt.ylabel("nonLinearDet")
    plt.legend()
  if order>=2 and p[-3]<0:
    log.warn("corrNonlinGetPar: consistency problem, second order coefficient should \
    be > 0, please double check result (plot=True) or try inverting the data and the\
    correct arguments")

  if returnCorrectedDet:
    return corrNonlin(nonLinearDet,p,data_0=data_0,correct_0=correct_0)
  else:
    return p

def plot_fit(deg, err):
    try:
        polcoefs = np.polyfit(x_rand, y_rand, deg)
    except np.RankWarning:
        pass
    if err:
        fig = plt.figure(figsize=(15, 10)) 
        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
        ax = plt.subplot(gs[0])
    else:
        fig, ax = plt.subplots(1, 1, figsize=(15, 10))
    # data
    ax.plot(x_cont, np.polyval(polcoefs, x_cont))
    ax.legend(['degree {}'.format(deg)])
    plot_dat(ax, err)
    # error
    if err:
        train_error.append(sum((np.polyval(polcoefs, x_rand)-y_rand)**2)/len(x_rand))
        test_error.append(sum((np.polyval(polcoefs, x_test)-y_test)**2)/len(x_test))
        ax = plt.subplot(gs[1])
        ax.plot(deg_list[0:len(train_error)], train_error, marker='o', c=(0, 0, 1))
        ax.plot(deg_list[0:len(test_error)], test_error, marker='o', c=(0.9, 0, 0))
        ax.set_xlabel('degree')
        ax.set_ylabel('error')
        ax.legend(['train', 'test'])
        ax.set_xlim(0.9, deg_list[-1]+0.1)
        ax.set_ylim(0, None)
        ax.set_yticks([])
        # layout
        plt.tight_layout()

    def update(self, rho):
        """
        Calculate the probability function for the given state of an harmonic
        oscillator (as density matrix)
        """

        if isket(rho):
            rho = ket2dm(rho)

        self.data = np.zeros(len(self.xvecs[0]), dtype=complex)
        M, N = rho.shape

        for m in range(M):
            k_m = pow(self.omega / pi, 0.25) / \
                sqrt(2 ** m * factorial(m)) * \
                exp(-self.xvecs[0] ** 2 / 2.0) * \
                np.polyval(hermite(m), self.xvecs[0])

            for n in range(N):
                k_n = pow(self.omega / pi, 0.25) / \
                    sqrt(2 ** n * factorial(n)) * \
                    exp(-self.xvecs[0] ** 2 / 2.0) * \
                    np.polyval(hermite(n), self.xvecs[0])

                self.data += np.conjugate(k_n) * k_m * rho.data[m, n]

def rssmodelwave(grating,grang,artic,cbin,cols):
#   compute wavelengths from model (this can probably be done using pyraf spectrograph model)
    spec=np.loadtxt(datadir+"spec.txt",usecols=(1,))
    Grat0,Home0,ArtErr,T2Con,T3Con=spec[0:5]
    FCampoly=spec[5:11]
    grname=np.loadtxt(datadir+"gratings.txt",dtype=str,usecols=(0,))
    grlmm,grgam0=np.loadtxt(datadir+"gratings.txt",usecols=(1,2),unpack=True)

    grnum = np.where(grname==grating)[0][0]
    lmm = grlmm[grnum]
    alpha_r = np.radians(grang+Grat0)
    beta0_r = np.radians(artic*(1+ArtErr)+Home0)-alpha_r
    gam0_r = np.radians(grgam0[grnum])
    lam0 = 1e7*np.cos(gam0_r)*(np.sin(alpha_r) + np.sin(beta0_r))/lmm
    ww = lam0/1000. - 4.
    fcam = np.polyval(FCampoly,ww)
    disp = (1e7*np.cos(gam0_r)*np.cos(beta0_r)/lmm)/(fcam/.015)
    dfcam = 3.162*disp*np.polyval([FCampoly[x]*(5-x) for x in range(5)],ww)
    T2 = -0.25*(1e7*np.cos(gam0_r)*np.sin(beta0_r)/lmm)/(fcam/47.43)**2 + T2Con*disp*dfcam
    T3 = (-1./24.)*3162.*disp/(fcam/47.43)**2 + T3Con*disp
    T0 = lam0 + T2 
    T1 = 3162.*disp + 3*T3
    X = (np.array(range(cols))+1-cols/2)*cbin/3162.
    lam_X = T0+T1*X+T2*(2*X**2-1)+T3*(4*X**3-3*X)
    return lam_X

def get_cont(x,y,n=1,sl=1.,sh=5.):
	orilen = len(x)
	coef = np.polyfit(x,y,n)
	res = y - np.polyval(coef,x)
	IH = np.where(res>0)[0]
	IL = np.where(res<0)[0]
	dev = np.mean(res[IH])
	I = np.where((res>-sl*dev) & (res<sh*dev))[0]
	J1 = np.where(res<=-sl*dev)[0]
	J2 = np.where(res>=sh*dev)[0]
	J = np.unique(np.hstack((J1,J2)))
	cond = True
	if len(J)==0 or len(x)< .3*orilen:
		cond=False
	while cond:
		x = np.delete(x,J)
		y = np.delete(y,J)
		coef = np.polyfit(x,y,n)
		res = y - np.polyval(coef,x)
		IH = np.where(res>0)[0]
		IL = np.where(res<0)[0]
		dev = np.mean(res[IH])
		I = np.where((res>-sl*dev) & (res<sh*dev))[0]
		J1 = np.where(res<=-sl*dev)[0]
		J2 = np.where(res>=sh*dev)[0]
		J = np.unique(np.hstack((J1,J2)))
		cond = True
		if len(J)==0 or len(x)< .1*orilen:
			cond=False
	return coef

def get_ratio(sciw,rat,n=3):
	rat = scipy.signal.medfilt(rat,11)
	lori = len(sciw)
	coef = np.polyfit(sciw,rat,n)
	res = rat - np.polyval(coef,sciw)
	rms = np.sqrt(np.mean(res**2))
	I = np.where(res> 3*rms)[0]
	I2 = np.where(res< -3*rms)[0]
	I = np.sort(np.hstack((I,I2)))
	cond = True
	if len(I) == 0 or len(sciw) < .3 * lori:
		cond = False

	while cond:
		#imax = np.argmax(res**2)
		#sciw = np.delete(sciw,imax)
		#rat  = np.delete(rat,imax)
		sciw = np.delete(sciw,I)
		rat  = np.delete(rat,I)
		coef = np.polyfit(sciw,rat,n)
		res = rat - np.polyval(coef,sciw)
		rms = np.sqrt(np.mean(res**2))
		I = np.where(res> 3*rms)[0]
		I2 = np.where(res< -3*rms)[0]
		I = np.sort(np.hstack((I,I2)))
		if len(I) == 0 or len(sciw) < .3 * lori:
			cond = False

	return coef

def get_rats(ZO,ZI,ZF,pars):
	ords = []
	for i in range(sc.shape[1]):
		J1 = np.where(mw > sc[0,i,-1])[0]
		J2 = np.where(mw < sc[0,i,0])[0]
		if len(J1)>0 and len(J2)>0:
			ords.append(i)
	ords = np.array(ords)
	mf = get_full_model(pars[0],pars[1],pars[2],pars[3],RES_POW)
	tmodf = np.zeros((sc.shape[1],sc.shape[2]))
	tscif = np.zeros((sc.shape[1],sc.shape[2]))
	test_plot = np.zeros((4,sc.shape[1],sc.shape[2]))
	for i in ords:
		I = np.where((mw>sc[0,i,0]) & (mw<sc[0,i,-1]))[0]
		modw = mw[I]
		modf = mf[I]
		sciw = sc[0,i]
		scif = sc[3,i]/np.median(sc[3,i])
		modf = pixelization(modw,modf,sciw)
		#IMB = np.where(mask_bin[i]!=0)[0]
		#modf /= modf[IMB].mean()
		mscif = scipy.signal.medfilt(scif,11)
		rat = modf/mscif
		INF = np.where(mscif!=0)[0]
		coef = get_ratio(sciw[INF],rat[INF])
		scif = scif * np.polyval(coef,sciw)
		mscif = mscif * np.polyval(coef,sciw)
		coef = get_cont(sciw,mscif)
		scif = scif / np.polyval(coef,sciw)
		#plot(sciw,scif)
		coef = get_cont(sciw,modf)
		modf = modf / np.polyval(coef,sciw)
		#plot(sciw,modf)	
		tmodf[i] = modf
		tscif[i] = scif
		test_plot[0,i] = sc[0,i]
		test_plot[1,i] = scif
		test_plot[2,i] = modf
		test_plot[3,i] = mask_bin[i]
	#show()
	#print vcdx
	hdu = pyfits.PrimaryHDU(test_plot)
	os.system('rm example.fits')
	hdu.writeto('example.fits')
	rat = tscif/tmodf

	nejx = np.arange(100)/100.
	ratsout = []

	for i in range(len(ZI)):
		ejy = rat[ZO[i],ZI[i]:ZF[i]]
		ejx = np.arange(len(ejy))/float(len(ejy))
		tck = interpolate.splrep(ejx,ejy,k=3)
		if len(ratsout)==0:
			ratsout = interpolate.splev(nejx,tck)
		else:
			ratsout = np.vstack((ratsout,interpolate.splev(nejx,tck)))
			#plot(interpolate.splev(nejx,tck))
	#show()
	return ratsout

def detrend(var, ax=None, lcopy=True, ldetrend=True, ltrend=False, degree=1, rcond=None, w=None,  
            lsmooth=False, lresidual=False, window_len=11, window='hanning'): 
  ''' subtract a linear trend from a time-series array (operation is in-place) '''
  # check input
  if not isinstance(var,np.ndarray): raise NotImplementedError # too many checks
  if lcopy: var = var.copy() # make copy - not in-place!
  # fit over entire array (usually not what we want...)
  if ax is None and ldetrend: ax = np.arange(var.size) # make dummy axis, if necessary
  if var.ndim != 1:
    shape = var.shape 
    var = var.ravel() # flatten array, if necessary
  else: shape = None
  # apply optional detrending
  if ldetrend or ltrend:
    # fit linear trend
    trend = np.polyfit(ax, var, deg=degree, rcond=rcond, w=w, full=False, cov=False)
    # evaluate and subtract linear trend
    if ldetrend and ltrend: raise ArgumentError("Can either return trend/polyfit or residuals, not both.")
    elif ldetrend and not ltrend: var -= np.polyval(trend, ax) # residuals
    elif ltrend and not ldetrend: var = np.polyval(trend, ax) # residuals
  # apply optional smoothing
  if lsmooth and lresidual: raise ArgumentError("Can either return smoothed array or residuals, not both.")
  elif lsmooth: var = smooth(var, window_len=window_len, window=window)  
  elif lresidual: var -= smooth(var, window_len=window_len, window=window)
  # return detrended and/or smoothed time-series
  if shape is not None: var = var.reshape(shape)
  return var

def dualPlot(age, meanWithin, meanBetween, title):

    fig, (within, between) = plt.subplots(1, 2, sharex=True, sharey=False)

    # fitshit
    wP = np.polyfit(age, meanWithin, 1)
    bP = np.polyfit(age, meanBetween, 1)
    xnew = np.arange(age.min() - 1, age.max() + 1, 0.1)
    wFit = np.polyval(wP, xnew)
    bFit = np.polyval(bP, xnew)

    within.set_title("within network")
    between.set_title("between network")

    withinCorr, withinP = st.pearsonr(age, meanWithin)
    within.plot(age, meanWithin, "k.")
    within.plot(xnew, wFit, "r", label=(str(np.round(withinCorr, 2)) + " " + str(np.round(withinP, 4))))
    within.set_xlabel("mean connectivity")
    within.set_ylabel("age")
    within.legend()

    betweenCorr, betweenP = st.pearsonr(age, meanBetween)
    between.plot(age, meanBetween, "k.")
    between.plot(xnew, bFit, "b", label=(str(np.round(betweenCorr, 2)) + " " + str(np.round(betweenP, 4))))
    between.set_xlabel("mean connectivity")
    between.set_ylabel("age")
    between.legend()

    fig.suptitle(title)
    plt.show()
    raw_input("Press Enter to continue...")
    plt.close()

def tfit(line):
    """
    Correct for temperature systematics.  Fit a polynomial to (teff,abund)
    and require that the corrected solar value be 0.  We cut on vsini, 

    returns:
    (fitabund,fitpar,t,abund)
    fitabund - the temperature corrected abundance
    fitpar   - the parameters to the polynomial fit
    t        - the temperature array
    abund    - the non-temp-corrected abundances

    """
    deg  = 3 # fit with a 3rd degree polynomial
    #define abundance for the particular line we're looking at
    p = getelnum.Getelnum(line)
    elstr = p.elstr

    conn = sqlite3.connect(os.environ['STARSDB'])
    cur = conn.cursor()

    #pull in the abundances and the non-corrected abundances
    cmd = 'SELECT '+elstr+'_abund_nt,teff FROM mystars WHERE '+globcut(elstr)
    cur.execute(cmd)
    arr = np.array(cur.fetchall() ) 
    abund,t = arr[:,0],arr[:,1]
    abund = abund - p.abnd_sol

    #fit the points
    fitpar = np.polyfit(t,abund,deg)
    #subtract out the fit, while requiring that the solar value be 0.
    fitpar[deg] = fitpar[deg] - np.polyval(fitpar,p.teff_sol)
    fitabund = abund - np.polyval(fitpar,t)
    return (fitabund,fitpar,t,abund)