Python numpy.polyval() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def plot(self, dataset, path, show=False):
        with PdfPages(path) as pdf:
            x_vals = dataset.data['T'].tolist()
            y_vals = dataset.data[self.symbol].tolist()
            plt.plot(x_vals, y_vals, 'ro', alpha=0.4, markersize=4)

            x_vals2 = np.linspace(min(x_vals), max(x_vals), 80)
            fx = np.polyval(self._coeffs, x_vals2)
            plt.plot(x_vals2, fx, linewidth=0.3, label='')

            plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 4))
            plt.legend(loc=3, bbox_to_anchor=(0, 0.8))
            plt.title('$%s$ vs $T$' % self.display_symbol)
            plt.xlabel('$T$ (K)')

            plt.ylabel('$%s$ (%s)' % (self.display_symbol, self.units))

            fig = plt.gcf()
            pdf.savefig(fig)
            plt.close()

        if show:
            webbrowser.open_new(path) 

Example 2

def calibrate(self):
        # generate sequence
        self.set()
        # run and load normalized data
        data, _ = self.run(norm_pts = {self.qubit_names[0]: (0, 1), self.qubit_names[1]: (0, 2)})
        # select target qubit
        data_t = data[self.qubit_names[1]]
        # fit
        self.opt_par, all_params_0, all_params_1 = fit_CR([self.lengths, self.phases, self.amps], data_t, self.cal_type)
        # plot the result
        xaxis = self.lengths if self.cal_type==CR_cal_type.LENGTH else self.phases if self.cal_type==CR_cal_type.PHASE else self.amps
        finer_xaxis = np.linspace(np.min(xaxis), np.max(xaxis), 4*len(xaxis))
        self.plot["Data 0"] = (xaxis,       data_t[:len(data_t)//2])
        self.plot["Fit 0"] =  (finer_xaxis, np.polyval(all_params_0, finer_xaxis) if self.cal_type == CR_cal_type.AMP else sinf(finer_xaxis, *all_params_0))
        self.plot["Data 1"] = (xaxis,       data_t[len(data_t)//2:])
        self.plot["Fit 1"] =  (finer_xaxis, np.polyval(all_params_1, finer_xaxis) if self.cal_type == CR_cal_type.AMP else sinf(finer_xaxis, *all_params_1))
        return (str.lower(self.cal_type.name), self.opt_par) 

Example 3

def band_center(spectrum, low_endmember=None, high_endmember=None, degree=3):
    x = spectrum.index
    y = spectrum

    if not low_endmember:
        low_endmember = x[0]
    if not high_endmember:
        high_endmember = x[-1]

    ny = y[low_endmember:high_endmember]

    fit = np.polyfit(ny.index, ny, degree)

    center_fit = Series(np.polyval(fit, ny.index), ny.index)
    center = band_minima(center_fit)

    return center, center_fit 

Example 4

def fit_cubic(y0, y1, g0, g1):
    """Fit cubic polynomial to function values and derivatives at x = 0, 1.

    Returns position and function value of minimum if fit succeeds. Fit does
    not succeeds if

    1. polynomial doesn't have extrema or
    2. maximum is from (0,1) or
    3. maximum is closer to 0.5 than minimum
    """
    a = 2*(y0-y1)+g0+g1
    b = -3*(y0-y1)-2*g0-g1
    p = np.array([a, b, g0, y0])
    r = np.roots(np.polyder(p))
    if not np.isreal(r).all():
        return None, None
    r = sorted(x.real for x in r)
    if p[0] > 0:
        maxim, minim = r
    else:
        minim, maxim = r
    if 0 < maxim < 1 and abs(minim-0.5) > abs(maxim-0.5):
        return None, None
    return minim, np.polyval(p, minim) 

Example 5

def __call__(self, dispersion, *parameters):
        """
        Generate data at the dispersion points, given the parameters.

        :param dispersion:
            An array of dispersion points to calculate the data for.

        :param parameters:
            Keyword arguments of the model parameters and their values.
        """

        function, profile_parameters = self._profiles[self.metadata["profile"]]

        N = len(profile_parameters)
        y = 1.0 - function(dispersion, *parameters[:N])

        # Assume rest of the parameters are continuum coefficients.
        if parameters[N:]:
            y *= np.polyval(parameters[N:], dispersion)
        
        return y 

Example 6

def meyeraux(x):
    """meyer wavelet auxiliary function.

    v(x) = 35*x^4 - 84*x^5 + 70*x^6 - 20*x^7.

    Parameters
    ----------
    x : array
        grid points

    Returns
    -------
    y : array
        values at x
    """
    # Auxiliary def values.
    y = np.polyval([-20, 70, -84, 35, 0, 0, 0, 0], x) * (x >= 0) * (x <= 1)
    y += (x > 1)
    return y 

Example 7

def get_bias_coeff_from_T(self, master_bias_temp, master_bias_level,
                              frame_temp, calibrated_params):
        """
        :param master_bias_temp: Temperature of the master bias frame.

        :param master_bias_level: Median of the master bias frame.

        :param frame_temp: Temperature of the frame to correct.

        :param calibrated_params: parameters [a, b] of the
          function bias_level(T) = aT + b. T is in degrees C and
          bias_level(T) is the median of the bias frame at the
          given temperature.
        """
        calib_master_bias_level = np.polyval(
            calibrated_params, master_bias_temp)
        calib_frame_bias_level = np.polyval(
            calibrated_params, frame_temp)
        return  calib_frame_bias_level / calib_master_bias_level 

Example 8

def csalbr(tau):
    # Previously 3 functions csalbr fintexp1, fintexp3
    a = [-.57721566, 0.99999193, -0.24991055, 0.05519968, -0.00976004,
         0.00107857]
    #xx = a[0] + a[1] * tau + a[2] * tau**2 + a[3] * tau**3 + a[4] * tau**4 + a[5] * tau**5
    #xx = np.polyval(a[::-1], tau)

    # xx = a[0]
    # xftau = 1.0
    # for i in xrange(5):
    #     xftau = xftau*tau
    #     xx = xx + a[i] * xftau
    fintexp1 = np.polyval(a[::-1], tau) - np.log(tau)
    fintexp3 = (np.exp(-tau) * (1.0 - tau) + tau**2 * fintexp1) / 2.0

    return (3.0 * tau - fintexp3 *
            (4.0 + 2.0 * tau) + 2.0 * np.exp(-tau)) / (4.0 + 3.0 * tau)

# From crefl.1.7.1 

Example 9

def fit_dimensions(self, data, fit_TTmax = True):
        """ 
        if fit_max is True, use blade length profile to adjust dHS_max
        """
        
        if (fit_TTmax and 'L_blade' in data.columns):
            dat = data.dropna(subset=['L_blade'])
            xn = self.xn(dat['rank'])
            fit = numpy.polyfit(xn,dat['L_blade'],7)
            x = numpy.linspace(min(xn), max(xn), 500)
            y = numpy.polyval(fit,x)
            self.xnmax = x[numpy.argmax(y)]
            self.TTmax = self.TTxn(self.xnmax)
            

        self.scale = {k: numpy.mean(data[k] / self.predict(k,data['rank'], data['nff'])) for k in data.columns if k in self.ref.columns}

        return self.scale 

Example 10

def __init__(self, points, values, *args, **kwds):
        """
        Parameters
        ----------
        points : nd array (npoints, ndim)
        values : 1d array (npoints,)
        **kwds : keywords to [avg]polyfit()
        """
        self.points = self._fix_shape_init(points)
        assert self.points.ndim == 2, "points is not 2d array"
        self.values = values
        if self._has_keys(kwds, ['degrange', 'degmin', 'degmax', 'levels']):
            self.fitfunc = avgpolyfit
            self.evalfunc = avgpolyval
        else:
            self.fitfunc = polyfit
            self.evalfunc = polyval
        self.fit = self.fitfunc(self.points, self.values, *args, **kwds) 

Example 11

def vp2dp(vp,p=[],temp=[],enhance=False):
    """
    Convert a volume mixing ratio to a dew point ( and vapour pressure )
    Using ITS-90 correction of Wexler's formula
    Optional enhancement factors for non ideal 
    """  
    vp=np.atleast_1d(vp)
    c=np.array([-9.2288067e-06, 0.46778925, -20.156028, 207.98233],dtype='f8')
    d=np.array([-7.5172865e-05, 0.0056577518, -0.13319669, 1],dtype='f8')
    lnes=np.log(vp*1e2)
    dp=np.polyval(c,lnes)/np.polyval(d,lnes)
    if(enhance and len(p)>0) :
        if(len(temp)==0):
            temp=dp
        A=np.array([8.5813609e-09, -6.7703064e-06, 0.001807157, -0.16302041],dtype='f8')
        B=np.array([6.3405286e-07, -0.00077326396, 0.34378043, -59.890467],dtype='f8')
        alpha=np.polyval(A,temp)
        beta=np.exp(np.polyval(B,temp))
        ef=np.exp(alpha*(1-vp/p)+beta*(p/vp-1))
        vp=vp/ef
        lnes=np.log(vp*1e2)
        dp=np.polyval(c,lnes)/np.polyval(d,lnes)
    return dp 

Example 12

def eval_trace_poly(self, use_poly=False, smoothing_length=25):
        sel = self.trace_x != self.flag
        if use_poly:
            self.trace = np.polyval(self.trace_polyvals, 
                                    1. * np.arange(self.D) / self.D)
        else:
            self.trace = np.zeros((self.D,))
            init_x = np.where(sel)[0][0]
            fin_x = np.where(sel)[0][-1]
            self.trace[init_x:fin_x] = np.interp(np.arange(init_x,fin_x), 
                                                self.trace_x[sel],
                                                self.trace_y[sel])
            self.trace[init_x:fin_x] = biweight_filter(self.trace[init_x:fin_x], 
                                              smoothing_length)
            ix = int(init_x+smoothing_length/2+1)
            fx = int(init_x+smoothing_length/2+1 + smoothing_length*2)
            p1 = np.polyfit(np.arange(ix,fx), self.trace[ix:fx], 1)
            self.trace[:ix] = np.polyval(p1, np.arange(ix))
            ix = int(fin_x-smoothing_length/2-1 - smoothing_length*2)
            fx = int(fin_x-smoothing_length/2) 
            pf = np.polyfit(np.arange(ix,fx), self.trace[ix:fx], 1)
            self.trace[fx:self.D] = np.polyval(pf, np.arange(fx,self.D)) 

Example 13

def calculate_wavelength_new(x, solar_spec, fibers, fibn, group,
                              smooth_length=21, init_lims=None, order=3, 
                              init_sol=None, debug=False, interactive=False, 
                              nbins=21, wavebuff=100, plotbuff=85, 
                              fixscale=True, use_leastsq=False, res=1.9):
    L = len(x)
    if init_lims is None:
        init_lims = [np.min(solar_spec[:,0]), np.max(solar_spec[:,0])]
    if init_sol is not None:
        init_wave_sol = np.polyval(init_sol, 1. * x / L)
    y_sun = solar_spec[:,1]  
    lowfib = np.max([0,fibn-group])
    highfib = np.min([len(fibers)-1,fibn+group])
    y = np.array([biweight_filter(fibers[i].spectrum, smooth_length)
                  / fibers[i].spectrum 
                  for i in xrange(lowfib,highfib)])
    y = biweight_location(y,axis=(0,))
    bins = np.linspace(init_lims[0], init_lims[1], nbins)
    bins = bins[1:-1]
    scale = 1.*(init_lims[1] - init_lims[0])/L
    wv0 = init_lims[0]
    wave0_save = [] 
    scale_save = []
    x_save = []
    wave_save = [] 

Example 14

def calc_omega(cp):
    cp.insert
    a=[]
    for i in range(len(cp)):
        ptmp = []
        tmp = 0
        for j in range(len(cp)):
            if j != i:
                row = []
                row.insert(0,1/(cp[i]-cp[j]))
                row.insert(1,-cp[j]/(cp[i]-cp[j]))
                ptmp.insert(tmp,row)
                tmp += 1
        p=[1]
        for j in range(len(cp)-1):
            p = conv(p,ptmp[j])
        pint = numpy.polyint(p)
        arow = []
        for j in range(len(cp)):
            arow.append(numpy.polyval(pint,cp[j]))
        a.append(arow)
    return a 

Example 15

def calc_adot(cp,order=1):
    a=[]
    for i in range(len(cp)):
        ptmp = []
        tmp = 0
        for j in range(len(cp)):
            if j != i:
                row = []
                row.insert(0,1/(cp[i]-cp[j]))
                row.insert(1,-cp[j]/(cp[i]-cp[j]))
                ptmp.insert(tmp,row)
                tmp += 1
        p=[1]
        for j in range(len(cp)-1):
            p = conv(p,ptmp[j])
        pder = numpy.polyder(p,order)
        arow = []
        for j in range(len(cp)):
            arow.append(numpy.polyval(pder,cp[j]))
        a.append(arow)
    return a 

Example 16

def calc_afinal(cp):
    afinal=[]
    for i in range(len(cp)):
        ptmp = []
        tmp = 0
        for j in range(len(cp)):
            if j != i:
                row = []
                row.insert(0,1/(cp[i]-cp[j]))
                row.insert(1,-cp[j]/(cp[i]-cp[j]))
                ptmp.insert(tmp,row)
                tmp += 1
        p=[1]
        for j in range(len(cp)-1):
            p = conv(p,ptmp[j])
        afinal.append(numpy.polyval(p,1.0))
    return afinal 

Example 17

def test_polyfit(self):
        c = np.array([3., 2., 1.])
        x = np.linspace(0, 2, 7)
        y = np.polyval(c, x)
        err = [1, -1, 1, -1, 1, -1, 1]
        weights = np.arange(8, 1, -1)**2/7.0

        # check 1D case
        m, cov = np.polyfit(x, y+err, 2, cov=True)
        est = [3.8571, 0.2857, 1.619]
        assert_almost_equal(est, m, decimal=4)
        val0 = [[2.9388, -5.8776, 1.6327],
                [-5.8776, 12.7347, -4.2449],
                [1.6327, -4.2449, 2.3220]]
        assert_almost_equal(val0, cov, decimal=4)

        m2, cov2 = np.polyfit(x, y+err, 2, w=weights, cov=True)
        assert_almost_equal([4.8927, -1.0177, 1.7768], m2, decimal=4)
        val = [[8.7929, -10.0103, 0.9756],
               [-10.0103, 13.6134, -1.8178],
               [0.9756, -1.8178, 0.6674]]
        assert_almost_equal(val, cov2, decimal=4)

        # check 2D (n,1) case
        y = y[:, np.newaxis]
        c = c[:, np.newaxis]
        assert_almost_equal(c, np.polyfit(x, y, 2))
        # check 2D (n,2) case
        yy = np.concatenate((y, y), axis=1)
        cc = np.concatenate((c, c), axis=1)
        assert_almost_equal(cc, np.polyfit(x, yy, 2))

        m, cov = np.polyfit(x, yy + np.array(err)[:, np.newaxis], 2, cov=True)
        assert_almost_equal(est, m[:, 0], decimal=4)
        assert_almost_equal(est, m[:, 1], decimal=4)
        assert_almost_equal(val0, cov[:, :, 0], decimal=4)
        assert_almost_equal(val0, cov[:, :, 1], decimal=4) 

Example 18

def var(x):
	"""
	Creates a polynomial that describes the kernel width
	"""
	p=[0.2,0.02]
	return np.polyval(p,x)

#x sampling: 

Example 19

def extrap1d_constrained_linear_regression(x, y, xEval, side='right', numPts=10):
    """Perform extrapolation using constrained linear regression on part of the data (x,y).  Use numPts
    from either the left or right side of the data (specified by input variable side) as the input data.
    The linear regression is constrained to pass through the final point (x0, y0) (rightmost point if
    side=='right', leftmost if side=='left').  Data MUST be sorted.
    
    Inputs:
      x                 independent variable on the smaller domain (array)
      y                 dependent variable on the smaller domain (array)
      xEval             values of x at which to evaluate the linear regression model
      side              side of the data from which to perform linear regression ('left' or 'right')
      numPts            number of points to use in the linear regression (scalar)
    Outputs:
      yEval             values of dependent variable in the linear regression model evaluated at x_eval (array)
    """
    assert side=='left' or side=='right'
    if side=='left':
        xSide = x[:numPts]
        ySide = y[:numPts]
        x0 = x[0]
        y0 = y[0]
    elif side=='right':
        xSide = x[-numPts:]
        ySide = y[-numPts:]
        x0 = x[-1]
        y0 = y[-1]
        
    a = least_squares_slope(xSide, ySide, x0, y0)  # determine model (a, x0, y0)
    b = y0 - a*x0
    #y_eval = scipy.polyval([a,b], x_eval)
    yEval = a*(xEval - x0) + y0 # evaluate model on desired points
    return yEval 

Example 20

def remove_continuum_blockvisibility(vis: BlockVisibility, degree=1, mask=None) -> BlockVisibility:
    """ Fit and remove continuum visibility

    Fit a polynomial in frequency of the specified degree where mask is True
  
    :param vis:
    :param degree: Degree of polynomial
    :param mask:
    :return:
    """
    assert isinstance(vis, Visibility) or isinstance(vis, BlockVisibility), vis
    
    if mask is not None:
        assert numpy.sum(mask) > 2 * degree, "Insufficient channels for fit"
    
    nchan = len(vis.frequency)
    x = (vis.frequency - vis.frequency[nchan // 2]) / (vis.frequency[0] - vis.frequency[nchan // 2])
    for row in range(vis.nvis):
        for ant2 in range(vis.nants):
            for ant1 in range(vis.nants):
                for pol in range(vis.polarisation_frame.npol):
                    wt = numpy.sqrt(vis.data['weight'][row, ant2, ant1, :, pol])
                    if mask is not None:
                        wt[mask] = 0.0
                    fit = numpy.polyfit(x, vis.data['vis'][row, ant2, ant1, :, pol], w=wt, deg=degree)
                    prediction = numpy.polyval(fit, x)
                    vis.data['vis'][row, ant2, ant1, :, pol] -= prediction
    return vis 

Example 21

def remove_continuum_image(im: Image, degree=1, mask=None):
    """ Fit and remove continuum visibility in place
    
    Fit a polynomial in frequency of the specified degree where mask is True

    :param im:
    :param deg:
    :param mask:
    :return:
    """
    assert isinstance(im, Image)

    if mask is not None:
        assert numpy.sum(mask) > 2 * degree, "Insufficient channels for fit"

    nchan, npol, ny, nx = im.shape
    channels = numpy.arange(nchan)
    with warnings.catch_warnings():
        warnings.simplefilter('ignore')
        frequency = im.wcs.sub(['spectral']).wcs_pix2world(channels, 0)[0]
    frequency -= frequency[nchan // 2]
    frequency /= numpy.max(frequency)
    wt = numpy.ones_like(frequency)
    if mask is not None:
        wt[mask] = 0.0

    for pol in range(npol):
        for y in range(ny):
            for x in range(nx):
                    fit = numpy.polyfit(frequency, im.data[:, pol, y, x], w=wt, deg=degree)
                    prediction = numpy.polyval(fit, frequency)
                    im.data[:, pol, y, x] -= prediction
    return im 

Example 22

def phaserate(plotar, ms2mappings):
    if plotar.plotType!='phatime':
        raise RuntimeError, "phaserate() cannot run on plot type {0}".format( plotar.plotType )
    spm = ms2mappings.spectralMap
    # iterate over all plots and all data sets within the plots
    for k in plotar.keys():
        for d in plotar[k].keys():
            # get a reference to the data set
            dsref = plotar[k][d]
            # get the full data set label - we have access to all the data set's properties (FQ, SB, POL etc)
            n     = plots.join_label(k, d)
            # fit a line through the unwrapped phase 
            unw    = numpy.unwrap(numpy.deg2rad(dsref.yval))
            coeffs = numpy.polyfit(dsref.xval, unw, 1)
            # evaluate the fitted polynomial at the x-loci
            extray = numpy.polyval(coeffs, dsref.xval)
            # here we could compute the reliability of the fit
            diff   = unw - extray
            ss_tot = numpy.sum(numpy.square(unw - unw.mean()))
            ss_res = numpy.sum(numpy.square(diff))
            r_sq   = 1.0 - ss_res/ss_tot
            # compare std deviation and variance in the residuals after fit
            std_r  = numpy.std(diff)
            var_r  = numpy.var(diff)
            f      = spm.frequencyOfFREQ_SB(n.FQ, n.SB)
            rate   = coeffs[0]
            if var_r<std_r:
                print "{0}: {1:.8f} ps/s @ {2:5.4f}MHz [R2={3:.3f}]".format(n, rate/(2.0*numpy.pi*f*1.0e-12), f/1.0e6, r_sq )
                # before plotting wrap back to -pi,pi and transform to degrees
                dsref.extra = [ drawline_fn(n, dsref.xval, numpy.rad2deg(do_wrap(extray))) ] 

Example 23

def phasedbg(plotar, ms2mappings):
    global cache
    if plotar.plotType!='phatime':
        raise RuntimeError, "phasedbg() cannot run on plot type {0}".format( plotar.plotType )
    store = len(cache)==0
    # iterate over all plots and all data sets within the plots
    for k in plotar.keys():
        for d in plotar[k].keys():
            # get a reference to the data set
            dsref = plotar[k][d]
            # get the full data set label - we have access to all the data set's properties (FQ, SB, POL etc)
            n     = plots.join_label(k, d)
            # fit a line through the unwrapped phase 
            unw    = numpy.unwrap(numpy.deg2rad(dsref.yval))
            #coeffs = numpy.polyfit(dsref.xval, unw, 1)
            coeffs = numpy.polyfit(xrange(len(dsref.yval)), unw, 1)
            # evaluate the fitted polynomial at the x-loci
            extray = numpy.polyval(coeffs, dsref.xval)
            # here we could compute the reliability of the fit
            diff   = unw - extray
            # compare std deviation and variance in the residuals after fit
            std_r  = numpy.std(diff)
            var_r  = numpy.var(diff)
            coeffs = numpy.rad2deg(coeffs)
            if var_r<std_r:
                # decide what to do
                if store:
                    cache[n] = coeffs
                else:
                    # check if current key exists in cache; if so
                    # do differencing
                    otherCoeffs = cache.get(n, None)
                    if otherCoeffs is None:
                        print "{0}: not found in cache".format( n )
                    else:
                        delta = otherCoeffs - coeffs
                        print "{0.BL} {0.SB} {0.P}: dRate={1:5.4f} dOff={2:4.1f}".format(n, delta[0], delta[1]+360.0 if delta[1]<0.0 else delta[1])
                # before plotting wrap back to -pi,pi and transform to degrees
                #dsref.extra = [ drawline_fn(n, dsref.xval, numpy.rad2deg(do_wrap(extray))) ]
    if not store:
        cache = {} 

Example 24

def estimate_tau0(T0, atm):

	T0_fit = np.array([600,1000])
	tau0_fit = np.array([5,0.1])
	x = 1000.0/T0_fit
	y = 1.0*np.log10(tau0_fit)
	coeff_low = np.polyfit(x,y,1)

	T0_fit = np.array([1000,1600])
	tau0_fit = np.array([0.1,1e-4])
	x = 1000.0/T0_fit
	y = 1.0*np.log10(tau0_fit)
	coeff_high = np.polyfit(x,y,1)

	atm_fit = 1.0*np.array([1, 10, 20])
	mtp_fit = 1.0*np.array([1, 0.1, 0.02])
	x = np.sqrt(atm_fit)
	y = np.log10(mtp_fit)
	coeff_p = np.polyfit(x,y,1)

	if T0 < 1000:
		tau_1atm = (10 ** np.polyval(coeff_low, 1000.0/T0))
	else:
		tau_1atm = (10 ** np.polyval(coeff_high, 1000.0/T0))

	mtp = 10 ** np.polyval(coeff_p, np.sqrt(atm))

	tau = tau_1atm * mtp
	#print 'tau = '+str(tau)
	tau_level = 10 ** (np.floor(np.log10(tau)))
	#print 'tau_level = '+str(tau_level)
	tau_rounded = tau_level * np.round(tau/tau_level * 1.0)/1.0
	#print 'tau_rounded = '+str(tau_rounded)
	#print tau_level
	#print str(tau_rounded)

	#print tau_rounded
	tau_rounded = max(1e-2, tau_rounded)
	return float(str(tau_rounded)) 

Example 25

def compose_functions(cls,
                          x,
                          number_of_compositions=1,
                          functions=(np.sin, np.exp, np.square, np.polyval,
                                     np.tan, ),
                          clip_big_values=True,
                          clip_value=1e6):
        """Compose functions from an iterable of functions.

        This is a helper function to cover a more real life scenario of
        plottings.

        Arguments:
            x (numpy.array): array for which composed values will be computed.
            number_of_compositions (int): number of compositions of functions.
            functions (tuple): an iterable of functions.
            clip_big_values (bool): whether or not to limit function extremes.
            clip_value (float): limit values for function composition.

        Returns:
            y (numpy.array): array of composed functions

        """
        i = 0
        y = x
        while i < number_of_compositions:
            func = np.random.choice(functions)
            if func == np.polyval:
                n_coefs = np.random.randint(0, 10)
                coefs = np.random.randint(-50, 50, size=n_coefs)
                y = func(coefs, x)
            else:
                y = func(y)
            if clip_big_values:
                y = np.clip(y, -clip_value, clip_value)
            i += 1
        return y

    # Goes with 'fast' parameters by default. 

Example 26

def plot_fit(self):
        """Plot the training data in X array along with decision boundary
        """
        from matplotlib import pyplot as plt
        x1 = np.linspace(self.table.min(), self.table.max(), 100)
        #reverse self.theta as it requires coeffs from highest degree to constant term
        x2 = np.polyval(np.poly1d(self.theta[::-1]),x1)
        plt.plot(x1, x2, color='r', label='decision boundary');
        plt.scatter(self.X[:, 1], self.X[:, 2], s=40, c=self.y, cmap=plt.cm.Spectral)
        plt.legend()
        plt.show() 

Example 27

def ployfit( y, x=None, order = 20 ):
    '''
    fit data (one-d array) by a ploynominal function
    return the fitted one-d array
    '''
    if x is None:
        x = range(len(y))
    pol = np.polyfit(x, y, order)
    return np.polyval(pol, x) 

Example 28

def construct_lanes(data, pts_per_lane):
    poly_x, poly_y = data[:,:4], data[:,4:]
    lane_x = np.vstack(
        [np.polyval(poly_x[t], np.linspace(0, 50, pts_per_lane))
         for t in xrange(poly_x.shape[0])])
    lane_y = np.vstack(
        [np.polyval(poly_y[t], np.linspace(0, 50, pts_per_lane))
         for t in xrange(poly_y.shape[0])])
    lane = np.hstack([lane_x, lane_y])
    nnz = lane_x.sum(axis=1) > 0.
    return lane, nnz 

Example 29

def test_polyfit(self):
        c = np.array([3., 2., 1.])
        x = np.linspace(0, 2, 7)
        y = np.polyval(c, x)
        err = [1, -1, 1, -1, 1, -1, 1]
        weights = np.arange(8, 1, -1)**2/7.0

        # check 1D case
        m, cov = np.polyfit(x, y+err, 2, cov=True)
        est = [3.8571, 0.2857, 1.619]
        assert_almost_equal(est, m, decimal=4)
        val0 = [[2.9388, -5.8776, 1.6327],
                [-5.8776, 12.7347, -4.2449],
                [1.6327, -4.2449, 2.3220]]
        assert_almost_equal(val0, cov, decimal=4)

        m2, cov2 = np.polyfit(x, y+err, 2, w=weights, cov=True)
        assert_almost_equal([4.8927, -1.0177, 1.7768], m2, decimal=4)
        val = [[8.7929, -10.0103, 0.9756],
               [-10.0103, 13.6134, -1.8178],
               [0.9756, -1.8178, 0.6674]]
        assert_almost_equal(val, cov2, decimal=4)

        # check 2D (n,1) case
        y = y[:, np.newaxis]
        c = c[:, np.newaxis]
        assert_almost_equal(c, np.polyfit(x, y, 2))
        # check 2D (n,2) case
        yy = np.concatenate((y, y), axis=1)
        cc = np.concatenate((c, c), axis=1)
        assert_almost_equal(cc, np.polyfit(x, yy, 2))

        m, cov = np.polyfit(x, yy + np.array(err)[:, np.newaxis], 2, cov=True)
        assert_almost_equal(est, m[:, 0], decimal=4)
        assert_almost_equal(est, m[:, 1], decimal=4)
        assert_almost_equal(val0, cov[:, :, 0], decimal=4)
        assert_almost_equal(val0, cov[:, :, 1], decimal=4) 

Example 30

def calculate(self, **state):
        """
        Calculate the material physical property at the specified temperature
        in the units specified by the object's 'property_units' property.

        :param T: [K] temperature

        :returns: physical property value
        """
        super().calculate(**state)
        return np.polyval(self._coeffs, state['T']) 

Example 31

def fit_quartic(y0, y1, g0, g1):
    """Fit constrained quartic polynomial to function values and erivatives at x = 0,1.

    Returns position and function value of minimum or None if fit fails or has
    a maximum. Quartic polynomial is constrained such that it's 2nd derivative
    is zero at just one point. This ensures that it has just one local
    extremum.  No such or two such quartic polynomials always exist. From the
    two, the one with lower minimum is chosen.
    """
    def g(y0, y1, g0, g1, c):
        a = c+3*(y0-y1)+2*g0+g1
        b = -2*c-4*(y0-y1)-3*g0-g1
        return np.array([a, b, c, g0, y0])

    def quart_min(p):
        r = np.roots(np.polyder(p))
        is_real = np.isreal(r)
        if is_real.sum() == 1:
            minim = r[is_real][0].real
        else:
            minim = r[(r == max(-abs(r))) | r == -max(-abs(r))][0].real
        return minim, np.polyval(p, minim)

    D = -(g0+g1)**2-2*g0*g1+6*(y1-y0)*(g0+g1)-6*(y1-y0)**2  # discriminant of d^2y/dx^2=0
    if D < 1e-11:
        return None, None
    else:
        m = -5*g0-g1-6*y0+6*y1
        p1 = g(y0, y1, g0, g1, .5*(m+np.sqrt(2*D)))
        p2 = g(y0, y1, g0, g1, .5*(m-np.sqrt(2*D)))
        if p1[0] < 0 and p2[0] < 0:
            return None, None
        [minim1, minval1] = quart_min(p1)
        [minim2, minval2] = quart_min(p2)
        if minval1 < minval2:
            return minim1, minval1
        else:
            return minim2, minval2 

Example 32

def __call__(self, x, model_dispersion, model_normalized_flux, *parameters):

        # Marginalize over all models.

        assert len(parameters) == len(self.parameter_names)

        # Continuum.
        O = self.metadata["continuum_order"]
        continuum = 1.0 if 0 > O \
            else np.polyval(parameters[-(O + 1):][::-1], model_dispersion)

        y = model_normalized_flux * continuum

        # Smoothing?
        try:
            index = self.parameter_names.index("smoothing")
        except ValueError:
            None
        else:
            kernel = abs(parameters[index])
            y = ndimage.gaussian_filter1d(y, kernel, axis=-1)


        # Redshift?
        try:
            index = self.parameter_names.index("redshift")
        except ValueError:
            z = 0
        else:
            v = parameters[index]
            z = v/299792.458 # km/s

        return np.interp(x, model_dispersion * (1 + z), y) 

Example 33

def FIT(p, seq):
    
    k = len(p)
    predict = polyval(p, k+1)
    
    return predict 

Example 34

def project_xy(xy_coords, pvec):

    # get cubic polynomial coefficients given
    #
    #  f(0) = 0, f'(0) = alpha
    #  f(1) = 0, f'(1) = beta

    alpha, beta = tuple(pvec[CUBIC_IDX])

    poly = np.array([
        alpha + beta,
        -2*alpha - beta,
        alpha,
        0])

    xy_coords = xy_coords.reshape((-1, 2))
    z_coords = np.polyval(poly, xy_coords[:, 0])

    objpoints = np.hstack((xy_coords, z_coords.reshape((-1, 1))))

    image_points, _ = cv2.projectPoints(objpoints,
                                        pvec[RVEC_IDX],
                                        pvec[TVEC_IDX],
                                        K, np.zeros(5))

    return image_points 

Example 35

def calcmlmags(self, lightcurve):
		"""
		Returns a "lc.mags"-like array made using the ml-parameters.
		It has the same size as lc.mags, and contains the microlensing to be added to them.
		The lightcurve object is not changed !
		
		For normal use, call getmags() from the lightcurve.
		
		Idea : think about only returning the seasons mags to speed it up ? Not sure if reasonable, as no seasons defined outside ?
		"""
		jds = lightcurve.jds[self.season.indices] # Is this already a copy ? It seems so. So no need for an explicit copy().
		# We do not need to apply shifts (i.e. getjds()), as anyway we "center" the jds.
		
		
		# Old method :
		if self.mltype == "poly":
		
			refjd = np.mean(jds)
			jds -= refjd # This is apparently safe, it does not shifts the lightcurves jds.
		
			allmags = np.zeros(len(lightcurve.jds))
			allmags[self.season.indices] = np.polyval(self.params, jds) # probably faster then +=
			return allmags
		
		
		# Legendre polynomials :
		if self.mltype == "leg":
		
			rjd = (np.max(jds) - np.min(jds))/2.0
			cjd = (np.max(jds) + np.min(jds))/2.0
			jds = (jds - cjd)/rjd
		
			allmags = np.zeros(len(lightcurve.jds))
		
			for (n, p) in enumerate(self.params):
				allmags[self.season.indices] += p * ss.legendre(n)(jds)
		
			return allmags 

Example 36

def test_polyfit(self):
        c = np.array([3., 2., 1.])
        x = np.linspace(0, 2, 7)
        y = np.polyval(c, x)
        err = [1, -1, 1, -1, 1, -1, 1]
        weights = np.arange(8, 1, -1)**2/7.0

        # check 1D case
        m, cov = np.polyfit(x, y+err, 2, cov=True)
        est = [3.8571, 0.2857, 1.619]
        assert_almost_equal(est, m, decimal=4)
        val0 = [[2.9388, -5.8776, 1.6327],
                [-5.8776, 12.7347, -4.2449],
                [1.6327, -4.2449, 2.3220]]
        assert_almost_equal(val0, cov, decimal=4)

        m2, cov2 = np.polyfit(x, y+err, 2, w=weights, cov=True)
        assert_almost_equal([4.8927, -1.0177, 1.7768], m2, decimal=4)
        val = [[8.7929, -10.0103, 0.9756],
               [-10.0103, 13.6134, -1.8178],
               [0.9756, -1.8178, 0.6674]]
        assert_almost_equal(val, cov2, decimal=4)

        # check 2D (n,1) case
        y = y[:, np.newaxis]
        c = c[:, np.newaxis]
        assert_almost_equal(c, np.polyfit(x, y, 2))
        # check 2D (n,2) case
        yy = np.concatenate((y, y), axis=1)
        cc = np.concatenate((c, c), axis=1)
        assert_almost_equal(cc, np.polyfit(x, yy, 2))

        m, cov = np.polyfit(x, yy + np.array(err)[:, np.newaxis], 2, cov=True)
        assert_almost_equal(est, m[:, 0], decimal=4)
        assert_almost_equal(est, m[:, 1], decimal=4)
        assert_almost_equal(val0, cov[:, :, 0], decimal=4)
        assert_almost_equal(val0, cov[:, :, 1], decimal=4) 

Example 37

def _ir_calibrate(self, data):
        """IR calibration
        """

        cwl = self._header['block5']["central_wave_length"][0] * 1e-6
        c__ = self._header['calibration']["speed_of_light"][0]
        h__ = self._header['calibration']["planck_constant"][0]
        k__ = self._header['calibration']["boltzmann_constant"][0]
        a__ = (h__ * c__) / (k__ * cwl)

        #b__ = ((2 * h__ * c__ ** 2) / (1.0e6 * cwl ** 5 * data.data)) + 1

        data.data[:] *= 1.0e6 * cwl ** 5
        data.data[:] **= -1
        data.data[:] *= (2 * h__ * c__ ** 2)
        data.data[:] += 1

        #Te_ = a__ / np.log(b__)

        data.data[:] = a__ / np.log(data.data)

        c0_ = self._header['calibration']["c0_rad2tb_conversion"][0]
        c1_ = self._header['calibration']["c1_rad2tb_conversion"][0]
        c2_ = self._header['calibration']["c2_rad2tb_conversion"][0]

        #data.data[:] = c0_ + c1_ * Te_ + c2_ * Te_ ** 2

        data.data[:] = np.polyval([c2_, c1_, c0_], data.data)

        data.mask[data.data < 0] = True
        data.mask[np.isnan(data.data)] = True 

Example 38

def _discretize(self, observation):
        if not self._is_discrete:
            buckets = np.array(observation, dtype=np.int)
            for i in range(len(observation)):
                buckets[i] = np.digitize(observation[i], self._separators[i])
            observation = np.polyval(buckets, self._bins)

        return observation 

Example 39

def test_polyfit(self):
        c = np.array([3., 2., 1.])
        x = np.linspace(0, 2, 7)
        y = np.polyval(c, x)
        err = [1, -1, 1, -1, 1, -1, 1]
        weights = np.arange(8, 1, -1)**2/7.0

        # check 1D case
        m, cov = np.polyfit(x, y+err, 2, cov=True)
        est = [3.8571, 0.2857, 1.619]
        assert_almost_equal(est, m, decimal=4)
        val0 = [[2.9388, -5.8776, 1.6327],
                [-5.8776, 12.7347, -4.2449],
                [1.6327, -4.2449, 2.3220]]
        assert_almost_equal(val0, cov, decimal=4)

        m2, cov2 = np.polyfit(x, y+err, 2, w=weights, cov=True)
        assert_almost_equal([4.8927, -1.0177, 1.7768], m2, decimal=4)
        val = [[8.7929, -10.0103, 0.9756],
               [-10.0103, 13.6134, -1.8178],
               [0.9756, -1.8178, 0.6674]]
        assert_almost_equal(val, cov2, decimal=4)

        # check 2D (n,1) case
        y = y[:, np.newaxis]
        c = c[:, np.newaxis]
        assert_almost_equal(c, np.polyfit(x, y, 2))
        # check 2D (n,2) case
        yy = np.concatenate((y, y), axis=1)
        cc = np.concatenate((c, c), axis=1)
        assert_almost_equal(cc, np.polyfit(x, yy, 2))

        m, cov = np.polyfit(x, yy + np.array(err)[:, np.newaxis], 2, cov=True)
        assert_almost_equal(est, m[:, 0], decimal=4)
        assert_almost_equal(est, m[:, 1], decimal=4)
        assert_almost_equal(val0, cov[:, :, 0], decimal=4)
        assert_almost_equal(val0, cov[:, :, 1], decimal=4) 

Example 40

def _calc_MC(self):
        """ Internal function for option valuation. See ``calc_px()`` for complete documentation.

        :Authors:
            Oleg Melnikov <[email protected]>
        """
        rng_seed, deg, n, m = self.px_spec.rng_seed, self.px_spec.deg, self.px_spec.nsteps, self.px_spec.npaths
        sp = self._LT_specs()
        dt, df = sp['dt'], sp['df_dt']

        S0, vol = self.ref.S0, self.ref.vol
        K, r, signCP = self.K, self.rf_r, self._signCP

        np.random.seed(rng_seed)
        norm_mtx = np.random.normal((r - 0.5 * vol ** 2) * dt, vol * math.sqrt(dt), (n + 1, m))
        S = S0 * np.exp(np.cumsum(norm_mtx, axis=0))
        S[0] = S0
        payout = np.maximum(signCP * (S - K), 0)
        v = np.copy(payout)  # terminal payouts

        # Least-Squares Monte Carlo (LSM):
        for i in range(n - 1, 0, -1):          # American Option Valuation by Backwards Induction
            rg = np.polyfit(S[i], v[i + 1] * df, deg)      # fit 5th degree polynomial to PV of current inner values
            C = np.polyval(rg, S[i])              # continuation values.
            v[i] = np.where(payout[i] > C, payout[i], v[i + 1] * df)  # exercise decision
        v[0] = v[1] * df

        v0 = np.mean(v[0])
        self.px_spec.add(px=v0, submethod='Least Squares Monte Carlo (LSM)')

        return self 

Example 41

def _calc_MC(self):
        """ Internal function for option valuation.

        :Authors:
            Yen-fei Chen <[email protected]>
        """
        _ = self.px_spec;   n, m, rng_seed, keep_hist, deg = _.nsteps, _.npaths, _.rng_seed, _.keep_hist, _.deg
        _ = self.ref;       S0, vol, q = _.S0, _.vol, _.q
        _ = self;           T, K, rf_r, net_r, sCP = _.T, _.K, _.rf_r, _.net_r, _.signCP
        _ = self._LT_specs(); u, d, p, df, dt = _['u'], _['d'], _['p'], _['df_dt'], _['dt']

        np.random.seed(rng_seed)

        # option_px = np.zeros((n + 1, m) ,'d')
        S = np.zeros((n + 1, m), 'd')  # stock price matrix
        S[0, :] = S0  # initial value

        # stock price paths
        for t in range(1, n+1):
            random = scipy.stats.norm.rvs(loc=0, scale=1, size=m)
            S[t, :] = S[t-1, :] * np.exp((rf_r - vol**2 / 2) * dt + vol * random * np.sqrt(dt))

        option_px = np.maximum(sCP*(S - K), 0)  # payoff when not shout
        final_payoff = np.repeat(S[-1, :], n+1, axis=0).reshape(m, n + 1)
        shout_px = np.maximum(sCP*(final_payoff.transpose() - K), sCP * (S - K))

        for t in range (n - 1, -1, -1):  # valuation process is similar to American option
            rg = np.polyfit(S[t, :], df * np.array(option_px[t + 1, :]), deg) # regression at time t
            C = np.polyval(rg, S[t, :])  # continuation values
            # exercise decision: shout or not shout
            option_px[t, :] = np.where(shout_px[t, :] > C, shout_px[t, :], option_px[t+1,:] * df)

        self.px_spec.add(px=float(np.mean(option_px[0, :])), sub_method='Hull p.609')
        return self 

Example 42

def op(i, j, x, y):
    p0, a, b, c = orthopy.line.recurrence_coefficients.jacobi(
            i, 0, 0,
            # standardization='monic'
            standardization='p(1)=(n+alpha over n)'
            )
    val1 = orthopy.line.tools.evaluate_orthogonal_polynomial(
            (x-y)/(x+y), p0, a, b, c
            )

    val1 = numpy.polyval(scipy.special.jacobi(i, 0, 0), (x-y)/(x+y))

    # treat x==0, y==0 separately
    if isinstance(val1, numpy.ndarray):
        idx = numpy.where(numpy.logical_and(x == 0, y == 0))[0]
        val1[idx] = numpy.polyval(scipy.special.jacobi(i, 0, 0), 0.0)
    else:
        if numpy.isnan(val1):
            val1 = numpy.polyval(scipy.special.jacobi(i, 0, 0), 0.0)

    p0, a, b, c = orthopy.line.recurrence_coefficients.jacobi(
            j, 2*i+1, 0,
            # standardization='monic'
            standardization='p(1)=(n+alpha over n)'
            )
    val2 = orthopy.line.tools.evaluate_orthogonal_polynomial(
            1-2*(x+y), p0, a, b, c
            )
    # val2 = numpy.polyval(scipy.special.jacobi(j, 2*i+1, 0), 1-2*(x+y))

    flt = numpy.vectorize(float)
    return flt(
        numpy.sqrt(2*i + 1) * val1 * (x+y)**i
        * numpy.sqrt(2*j + 2*i + 2) * val2
        ) 

Example 43

def test_eval(t, ref, tol=1.0e-14):
    n = 5
    p0, a, b, c = orthopy.line.recurrence_coefficients.legendre(
            n, 'monic', symbolic=True
            )
    value = orthopy.line.evaluate_orthogonal_polynomial(t, p0, a, b, c)

    assert value == ref

    # Evaluating the Legendre polynomial in this way is rather unstable, so
    # don't go too far with n.
    approx_ref = numpy.polyval(legendre(n, monic=True), t)
    assert abs(value - approx_ref) < tol
    return 

Example 44

def test_eval_vec(t, ref, tol=1.0e-14):
    n = 5
    p0, a, b, c = orthopy.line.recurrence_coefficients.legendre(
            n, 'monic', symbolic=True
            )
    value = orthopy.line.evaluate_orthogonal_polynomial(t, p0, a, b, c)

    assert (value == ref).all()

    # Evaluating the Legendre polynomial in this way is rather unstable, so
    # don't go too far with n.
    approx_ref = numpy.polyval(legendre(n, monic=True), t)
    assert (abs(value - approx_ref) < tol).all()
    return 

Example 45

def test_clenshaw(tol=1.0e-14):
    n = 5
    _, _, alpha, beta = \
        orthopy.line.recurrence_coefficients.legendre(n, 'monic')
    t = 1.0

    a = numpy.ones(n+1)
    value = orthopy.line.clenshaw(a, alpha, beta, t)

    ref = math.fsum([
            numpy.polyval(legendre(i, monic=True), t)
            for i in range(n+1)])

    assert abs(value - ref) < tol
    return 

Example 46

def polyval(fit, points, der=0, avg=False):
    """Evaluate polynomial generated by ``polyfit()`` on `points`.

    Parameters
    ----------
    fit, points : see polyfit()
    der : int, optional
        Derivative order. Only for 1D, uses np.polyder().
    avg : bool, optional
        Internal hack, only used by ``avgpolyval()``.
    
    Notes
    -----
    For 1D we provide "analytic" derivatives using np.polyder(). For ND, we
    didn't implement an equivalent machinery. For 2D, you might get away with
    fitting a bispline (see Interpol2D) and use it's derivs. For ND, try rbf.py's RBF
    interpolator which has at least 1st derivatives for arbitrary dimensions.

    See Also
    --------
    :class:`PolyFit`, :class:`PolyFit1D`, :func:`polyfit`
    """
    pscale, pmin = fit['pscale'], fit['pmin']
    vscale, vmin = fit['vscale'], fit['vmin']
    if der > 0:
        assert points.shape[1] == 1, "deriv only for 1d poly (ndim=1)"
        # ::-1 b/c numpy stores poly coeffs in reversed order
        dcoeffs = np.polyder(fit['coeffs'][::-1], m=der)
        return np.polyval(dcoeffs, (points[:,0] - pmin[0,0]) / pscale[0,0]) / \
            pscale[0,0]**der * vscale
    else:
        vand = vander((points - pmin) / pscale, fit['deg'])
        if avg:
            return np.dot(vand, fit['coeffs']) * vscale
        else:
            return np.dot(vand, fit['coeffs']) * vscale + vmin 

Example 47

def test_compare_numpy():
    x = np.sort(np.random.rand(10))
    y = np.random.rand(10)
    yy1 = np.polyval(np.polyfit(x, y, 3), x)
    for scale in [True,False]:
        yy2 = num.PolyFit1D(x, y, 3, scale=scale)(x)
        assert np.allclose(yy1, yy2) 

Example 48

def Legendre_matrix(N, ctheta):
    r"""Matrix of weighted Legendre Polynominals.

    Computes a matrix of weighted Legendre polynominals up to order N for
    the given angles

    .. math::

        L_n(\theta) = \frac{2n+1}{4 \pi} P_n(\theta)

    Parameters
    ----------
    N : int
        Maximum order.
    ctheta : (Q,) array_like
        Angles.

    Returns
    -------
    Lmn : ((N+1), Q) numpy.ndarray
        Matrix containing Legendre polynominals.
    """
    if ctheta.ndim == 0:
        M = 1
    else:
        M = len(ctheta)
    Lmn = np.zeros([N+1, M], dtype=complex)
    for n in range(N+1):
        Lmn[n, :] = (2*n+1)/(4*np.pi) * np.polyval(special.legendre(n), ctheta)

    return Lmn 

Example 49

def eval_fibmodel_poly(self, use_poly=False):
        self.fibmodel = np.zeros((self.D, len(self.binx)))
        for i in xrange(len(self.binx)):
            if use_poly:
                self.fibmodel[:,i] = np.polyval(self.fibmodel_polyvals[:,i],
                                                1.* np.arange(self.D) / self.D)
            else:
                self.fibmodel[:,i] = np.interp(np.arange(self.D), 
                                               self.fibmodel_x, 
                                               self.fibmodel_y[:,i]) 

Example 50

def eval_wave_poly(self):
        self.wavelength = np.polyval(self.wave_polyvals, 
                                     1.* np.arange(self.D) / self.D) 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注