Python numpy 模块，polyfit() 实例源码

def fit(self):
        # Quadratic fit to get an initial guess for the parameters. 
        # Thanks: https://github.com/materialsproject/pymatgen
        # -> pymatgen/io/abinitio/eos.py 
        a, b, c = np.polyfit(self.volume, self.energy, 2)
        v0 = -b/(2*a)
        e0 = a*v0**2 + b*v0 + c
        b0 = 2*a*v0
        b1 = 4  # b1 is usually a small number like 4
        if not self.volume.min() < v0 and v0 < self.volume.max():
            raise StandardError('The minimum volume of a fitted parabola is not in the input volumes')

        # need to use lst2dct and dct2lst here to keep the order of parameters 
        pp0_dct = dict(e0=e0, b0=b0, b1=b1, v0=v0)
        target = lambda pp, v: self.energy - self.func(v, self.func.lst2dct(pp))
        pp_opt, ierr = leastsq(target, 
                               self.func.dct2lst(pp0_dct), 
                               args=(self.volume,))
        self.params = self.func.lst2dct(pp_opt)

def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""

    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))

    # Initialize replicates: bs_slope_reps, bs_intercept_reps
    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        # noinspection PyTupleAssignmentBalance
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

def SLOPE(self, param):
        class Context:
            def __init__(self, N):
                self.N = N
                self.q = deque([], self.N)
                self.x = [i for i in range(self.N)]

            def handleInput(self, value):
                if len(self.q) < self.N:
                    self.q.append(value)
                    return 0
                self.q.append(value)
                z1 = np.polyfit(self.x, self.q, 1)
                return z1[0]
        ctx = Context(param[1])
        result = param[0].apply(ctx.handleInput)
        return result

def FORCAST(self, param):
        class Context:
            def __init__(self, N):
                self.N = N
                self.q = deque([], self.N)
                self.x = [i for i in range(self.N)]

            def handleInput(self, value):
                if len(self.q) < self.N:
                    self.q.append(value)
                    return np.NaN
                z1 = np.polyfit(self.x, self.q, 1)
                fn = np.poly1d(z1)
                y = fn(self.N + 1)
                self.q.append(value)
                return y

        ctx = Context(param[1])
        result = param[0].apply(ctx.handleInput)
        return result

#????

def minScalErr(stec,el,z,thisBias):
    """
    this determines the slope of the vTEC vs. Elevation line, which
    should be minimized in the minimum scalloping technique for
    receiver bias removal
    inputs:
        stec - time indexed Series of slant TEC values
        el - corresponding elevation values, also Series
        z - mapping function values to convert to vTEC from entire file, may
            contain nans, Series
        thisBias - the bias to be tested and minimized
    """

    intel=np.asarray(el[stec.index],int) # bin the elevation values into int
    sTEC=np.asarray(stec,float)
    zmap = z[stec.index]
    c=np.array([(i,np.average((sTEC[intel==i]-thisBias)
                              /zmap[intel==i])) for i in np.unique(intel) if i>30])

    return np.polyfit(c[:,0],c[:,1],1)[0]

def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""
    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))

    # Initialize replicates: bs_slope_reps, bs_intercept_reps
    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""

    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))

    # Initialize replicates: bs_slope reps, bs_intercept_reps
    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

def draw_bs_pairs(x, y, func, size=1):
    """Perform pairs bootstrap for linear regression."""

    # Set up array of indices to sample from: inds
    inds = np.arange(len(x))
    # Initialize replicates
    bs_replicates = np.empty(size)
    # bs_intercept_reps = ____

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_replicates[i] = func(bs_x, bs_y)
        # bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_replicates

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def calculateKDP(phiDP, ran):
    # smooth phiDP field and take derivative
    numRan = ran.shape[0]
    kdp = np.ma.masked_all(phiDP.shape)
    smoothPhiDP = smPhiDP(phiDP, ran)

    # take derivative of kdp field
    winLen = 11
    rprof = ran[0:winLen*2-1]/1000.

    for i in range(numRan-winLen*3):
        numvalid = smoothPhiDP[:,i:i+winLen*2-1].count(axis=1)
        max_numv = np.max(numvalid)
        if max_numv==(winLen*2-1):
            kdp[numvalid==(winLen*2-1),i+winLen] = 0.5*np.polyfit(rprof,
                smoothPhiDP[numvalid==(winLen*2-1),i:i+winLen*2-1].transpose(), 1)[0]
    return kdp

# function for creating color map
#----------------------------------

def hplot(name, h, err, xlab="", ylab="", hstr="h", rate=None, fit=True, fig=True):
    from matplotlib.pyplot import figure, loglog, xlabel, ylabel, legend, show
    if fig is True:
        figure()
    loglog(h, err, 's-', label=name)
    if fit:
        p = numpy.polyfit(numpy.log(h), numpy.log(err), 1)
        C = numpy.exp(p[1])
        alpha = p[0]
        alg = [C*n**alpha for n in h]
        loglog(h, alg, '--', label="%.3f*%s^(%.3f)" %(C, hstr, alpha))
    if rate and h[0] != 0:
        alg = [err[0]/(h[0]**rate)*n**rate for n in h]
        loglog(h, alg, '--', label="%s^(%.3f)" %(hstr, rate))
    xlabel(xlab)
    ylabel(ylab)
    legend(loc='best')

# --- create interpolation for different h and orders ---

def default1():
    numRadialPts = 144;
    rhoMin = 0.1
    rhoMax = 0.9
    rho = np.linspace(rhoMin, rhoMax, numRadialPts);
    qbar = 0.854 + 2.184 * rho**2;

    minorRadius = 0.28;  # a
    majorRadius = 0.71;  # R0
    r = rho * minorRadius;
    safetyFactor = qbar / np.sqrt(1 - (r/majorRadius)**2);

    # fit
    polynomialDegree = 3;
    p = np.polyfit(rho, safetyFactor, polynomialDegree)
    qCoeffs = p[::-1]
    return qCoeffs

def PCSCH3(self,row):
        actuals=[];read=[]
        resistance = float(self.tbl.item(row,0).text())
        for a in np.linspace(.2e-3,1.5e-3,20):
            actuals.append( self.I.DAC.setCurrent(a) )
            time.sleep(0.001)
            read.append (self.I.get_average_voltage('CH3',samples=5))
        read = np.array(read)/resistance 
        actuals = np.array(actuals)
        sread = read*1e3
        sactuals = actuals*1e3
        self.DacCurves['PCS'].setData(sactuals,sread-sactuals)
        self.tbl.item(row,1).setText(string.join(['%.3f'%a for a in sread-sactuals],' '))
        if np.any(abs(read-actuals)>10e-6):self.setSuccess(self.tbl.item(row,1),0)
        else: self.setSuccess(self.tbl.item(row,1),1)

        fitvals = np.polyfit(self.I.DAC.CHANS['PCS'].VToCode(read),self.I.DAC.CHANS['PCS'].VToCode(actuals),1)
        print (fitvals)
        if list(fitvals):
            self.PCS_SLOPE = fitvals[0] #slope
            self.PCS_OFFSET = fitvals[1]    #offset

def PCSCH3(self,row):
        actuals=[];read=[]
        resistance = float(self.tbl.item(row,0).text())
        for a in np.linspace(.2e-3,1.5e-3,20):
            actuals.append( self.I.DAC.setCurrent(a) )
            time.sleep(0.001)
            read.append (self.I.get_voltage('CH3',samples=5))
        read = np.array(read)/resistance 
        actuals = np.array(actuals)
        sread = read*1e3
        sactuals = actuals*1e3
        self.DacCurves['PCS'].setData(sactuals,sread-sactuals)
        self.tbl.item(row,1).setText(string.join(['%.3f'%a for a in sread-sactuals],' '))
        if np.any(abs(read-actuals)>10e-6):self.setSuccess(self.tbl.item(row,1),0)
        else: self.setSuccess(self.tbl.item(row,1),1)

        fitvals = np.polyfit(self.I.DAC.CHANS['PCS'].VToCode(read),self.I.DAC.CHANS['PCS'].VToCode(actuals),1)
        print (fitvals)
        if list(fitvals):
            self.PCS_SLOPE = fitvals[0] #slope
            self.PCS_OFFSET = fitvals[1]    #offset

def get_linear_regression(xdata, ydata):
    """ Calculate the coefficients of the linear equation corresponding
    to the linear regression of a series of points. """
    try:
        import numpy
        return tuple(numpy.polyfit(xdata, ydata, 1))
    except ImportError:
        # numpy not available
        # try something approximate and simple
        datasize = len(xdata)
        sum_x = float(sum(xdata))
        sum_y = float(sum(ydata))
        sum_xx = float(sum(map(lambda x: x * x, xdata)))
        sum_products = float(sum([xdata[i] * ydata[i]
                                  for i in range(datasize)]))
        a = (sum_products - sum_x * sum_y / datasize) / (
            sum_xx - (sum_x * sum_x) / datasize)
        b = (sum_y - a * sum_x) / datasize
        return a, b

def polyfit_baseline(bands, intensities, poly_order=5, num_stdv=3.,
                     max_iter=200):
    '''Iteratively fits a polynomial, discarding far away points as peaks.
    Similar in spirit to ALS and related methods.
    Automated method for subtraction of fluorescence from biological Raman spectra
    Lieber & Mahadevan-Jansen 2003
    '''
    fit_pts = intensities.copy()
    # precalculate [x^p, x^p-1, ..., x^1, x^0]
    poly_terms = bands[:, None] ** np.arange(poly_order, -1, -1)
    for _ in range(max_iter):
        coefs = np.polyfit(bands, fit_pts.T, poly_order)
        baseline = poly_terms.dot(coefs).T
        diff = fit_pts - baseline
        thresh = diff.std(axis=-1) * num_stdv
        mask = diff > np.array(thresh, copy=False)[..., None]
        unfitted = np.count_nonzero(mask)
        if unfitted == 0:
            break
        fit_pts[mask] = baseline[mask]  # these points are peaks, discard
    else:
        print("Warning: polyfit_baseline didn't converge in %d iters" % max_iter)
    return baseline

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

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def trendLine(self, axis_choose=None):
        stable_sec= int(self.record_sec_le.text())
        stable_count = int(stable_sec * (1/0.007))

        if axis_choose:
            axis = axis_choose
        else:
            axis = str(self.axis_combobox.currentText())

        x = self.raw_data['time'][:stable_count]
        y = self.raw_data[axis][:stable_count]
        coefficients = np.polyfit(x,y,1)
        p = np.poly1d(coefficients)
        coefficient_of_dermination = r2_score(y, p(x))

        self.trendLine_content1_label.setText("Trendline: " + str(p))
        self.trendLine_content2_label.setText("R: " + str(coefficient_of_dermination))
        return coefficients

def FitPolynomial(fig, x, y, n):

#solution, res = np.polynomial.polynomial.polyfit(x,y,order,full=True)

    solution, C_p = np.polyfit(x, y, n, cov=True)  # C_z is estimated covariance matrix

    # Do the interpolation for plotting:
    xrange = fig.get_xlim()
    t = np.linspace(xrange[0],xrange[1],100)
    # Matrix with rows 1, t, t**2, ...:
    TT = np.vstack([t**(n-i) for i in range(n+1)]).T
    yi = np.dot(TT, solution)  # matrix multiplication calculates the polynomial values
    C_yi = np.dot(TT, np.dot(C_p, TT.T)) # C_y = TT*C_z*TT.T
    sig_yi = np.sqrt(np.diag(C_yi))  # Standard deviations are sqrt of diagonal


    fig.plot(t,yi, 'k-')
    # fig.plot(t,yi+sig_yi, 'k--')
    # fig.plot(t,yi-sig_yi, 'k--')

    return solution

# Take the log of the y value

def compute_xvvr(self):
        """ Return xvv(r) matrix """
        r = np.array([i*self.dr for i in range(self.ngrid)])
        k = self.get_k()
        xvvr = [["" for i in range(self.nsites)] for j in range(self.nsites)]
        for i in range(self.nsites):
            for j in range(self.nsites):
                xvvk_ij = self.xvv_data[:,i,j]
                xvvr_ij = pubfft.sinfti(xvvk_ij*k, self.dr, -1)/r
#                n_pots_for_interp = 6
#                r_for_interp = r[1:n_pots_for_interp+1]
#                xvvr_for_interp = xvvr_ij[:n_pots_for_interp]
#                poly_coefs = np.polyfit(r_for_interp, xvvr_for_interp, 3)
#                poly_f = np.poly1d(poly_coefs)
#                xvvr[i][j] = [poly_f(0)]
                xvvr[i][j] = xvvr_ij
        return r, np.swapaxes(xvvr, 0, 2)

def compute_zr(self):
        """ Return z(r) matrix """
        r = np.array([i*self.dr for i in range(self.ngrid)])
        k, zk = self.compute_zk()
        print 'computed zk',zk.shape
        zr = [["" for i in range(self.nsites)] for j in range(self.nsites)]
        for i in range(self.nsites):
            for j in range(self.nsites):
                zk_ij = zk[1:,i,j]
                zr_ij = pubfft.sinfti(zk_ij*k[1:], self.dr, -1)/r[1:]
                #zr_ij = np.abs(fftpack.fft(zk_ij))
                n_pots_for_interp = 6
                r_for_interp = r[1:n_pots_for_interp+1]
                zr_for_interp = zr_ij[:n_pots_for_interp]
                poly_coefs = np.polyfit(r_for_interp, zr_for_interp, 3)
                poly_f = np.poly1d(poly_coefs)
                zr[i][j] = [poly_f(0)]
                zr[i][j].extend(zr_ij)
        return r, np.swapaxes(zr, 0, 2)

def smooth(x, y, weights):
    '''
    in case the NLF cannot be described by 
    a square root function
    commit bounded polynomial interpolation
    '''
    # Spline hard to smooth properly, therefore solfed with
    # bounded polynomal interpolation
    # ext=3: no extrapolation, but boundary value
#     return UnivariateSpline(x, y, w=weights,
#                             s=len(y)*weights.max()*100, ext=3)

#     return np.poly1d(np.polyfit(x,y,w=weights,deg=2))
    p = np.polyfit(x, y, w=weights, deg=2)
    if np.any(np.isnan(p)):
        # couldn't even do polynomial fit
        # as last option: assume constant noise
        my = np.average(y, weights=weights)
        return lambda x: my
    return lambda xint: np.poly1d(p)(np.clip(xint, x[0], x[-1]))

def __init__(self,xpixel,ypixel,library,zoomLevel):
        self.xpixel=xpixel
        self.ypixel=ypixel
        #load csv data. Rows are zoom levels, column are latitude
        #and values are meters per pixel
        with open(library,'rb') as csvfile:
            metersperpixel = list(csv.reader(csvfile,delimiter=","))
        latitudes=[]
        #convert to floats
        for element in metersperpixel[0][1:]:
            latitudes.append(float(element))
        for row in range(1,len(metersperpixel)):
            res_values=[]
            if int(metersperpixel[row][0])==zoomLevel:
            #convert to floats
                for element in metersperpixel[row][1:]:
                    res_values.append(float(element))
                self.fitvalues=\
                    numpy.polyfit(latitudes, #fit to latitutde values
                    res_values, 3)  #3rd degree polynomial fit
                print "Fit done."   
                break

def callback_click_ok(self,data):
        points = np.array(self.interal_data)    #[(1, 1), (2, 2), (3, 3), (7, 3), (9, 3)]
        # get x and y vectors
        x = points[:,0]
        y = points[:,1]

        # calculate polynomial
        terms=int(self.sp.value())
        self.ret = np.polyfit(x, y, terms)
        f = np.poly1d(self.ret)

        tot=""
        val=0
        for i in range(0,len(self.ret)):
            p=len(self.ret)-1-i
            tot=tot+str(self.ret[i])+"*pow(w,"+str(p)+")"+"+"
        tot=tot[:-1]
        self.ret_math=tot

        self.close()

def build_model_poly(detection_pairs, beacon_sdoa, nominal_sample_rate, deg=2):
    if len(detection_pairs) < deg + 1:
        # not enough beacon transmissions
        return None

    soa0 = np.array([d[0].soa for d in detection_pairs])
    soa1 = np.array([d[1].soa for d in detection_pairs])
    soa1at0 = soa1 + np.array(beacon_sdoa)
    coef = np.polyfit(soa1at0, soa0, deg)
    fit = np.poly1d(coef)
    # residuals = soa0 - fit(soa1at0)
    # print(np.mean(residuals))

    def evaluate(det0, det1):
        return (det0.soa - fit(det1.soa)) / nominal_sample_rate

    return evaluate

def ExpLine(*args, **kwargs):
        rate = kwargs.pop('rate', None)
        log_data = kwargs.pop('log_data', True)
        data = kwargs.pop('data', None)
        const = kwargs.pop('const', 1)
        if rate is None:
            assert data is not None and len(data) > 0, "rate or data must be given"
            x = np.array([d[0] for d in data])
            y = np.array([d[1] for d in data])
            if len(x) > 0 and len(y) > 0:
                if log_data:
                    rate = np.polyfit(np.log(x), np.log(y), 1)[0]
                else:
                    rate = np.polyfit(x, y, 1)[0]
        if "label" in kwargs:
            kwargs["label"] = kwargs.pop("label").format(rate=rate)
        return FunctionLine2D(*args, fn=lambda x, r=rate:
                              const*np.array(x)**r, data=data,
                              log_data=log_data, **kwargs)

def redo_fit(self):
        lx0, lx1 = self.power_plot_lr.getRegion()

        x0, x1 = 10**lx0, 10**lx1


        X = self.dat['power_meter_power']

        n = len(X)

        ii0 = np.argmin(np.abs(X[:n//2+1]-x0))
        ii1 = np.argmin(np.abs(X[:n//2+1]-x1))

        print(ii0,ii1)

        m, b = np.polyfit(np.log10(X[ii0:ii1]), np.log10(self.power_plot_y[ii0:ii1]), deg=1)
        print("fit", m,b) 

        fit_data = 10**(np.poly1d((m,b))(np.log10(X)))
        print("fit_data", fit_data)    
        self.power_fit_plotcurve.setData(X, fit_data)

        self.fit_text.setHtml("<h1>I<sup>{:1.2f}</sup></h1>".format(m))
        self.fit_text.setPos(0.5*(lx0+lx1), np.log10(fit_data[(ii0+ii1)//2]))

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def scatter_crimes_population():
    """creates a scatter plot using the values of Population and Crimes Per 100000.
    iterates through the database and reads all the data in the given headers and creates plots for each data point
    """
    x = df["Number of Crimes"].values
    y = df["Population"].values
    assert len(x) == len(y)
    df["Crimes Per 100000"] = np.array([(x[i] / y[i]) * 100000.0 for i in range(len(x))], dtype="float32")

    ax = df.plot.scatter(y="Population", x="Crimes Per 100000")

    for index, row in df.iterrows():
        ax.annotate(row["Community Name"], (row["Crimes Per 100000"], row["Population"]),
                    size=7,
                    color='darkslategrey')

    x = df["Crimes Per 100000"].values
    y = df["Population"].values

    m, b = np.polyfit(x, y, 1)
    line = plt.plot(x, m * x + b, 'b--')
    plt.setp(line, color='orange', alpha=0.5, linewidth=2.0)
    plt.show()

def linbleeched_F0(data):
    """
    Calculate a linear fit (:math:`y(t)=m t+y_0)` for each pixel, which is assumed to correct for bleeching effects.

    :param data: he video data of shape (M,N,T).
    :return: tuple (m,y0) with two images each with shape (M,N).
    """

    # generate c coordinates
    x = numpy.arange(data.shape[-1])
    # reshape the data to two d array, first dimension is pixel index, second dimension is time
    d = numpy.reshape(data, (data.shape[0] * data.shape[1], data.shape[-1]))
    # find fit parameters
    m, y0 = numpy.polyfit(x, d.T, 1)
    # reshape fit parameters back to image shape
    return m.reshape(data.shape[0:2]), y0.reshape(data.shape[0:2])

def RotatePoints(self, x_in, y_in):
        #Rotational matrix
        angle = np.pi/1.065
        rot = np.array([[np.cos(angle),-np.sin(angle)],[np.sin(angle),np.cos(angle)]])

        #Rotation Stuff
        x = 0.
        y = np.log(1000.)*-self.rd
        d1 = np.array([[x],[y]])
        V1 = np.dot(rot,d1)
        B1 = np.linspace(d1[0],V1[0],100)
        B2 = np.linspace(d1[1],V1[1],100)
        p = np.polyfit(B1,B2,1)

        #Skew-T y-axis
        y_out = -self.rd*np.log(y_in) 

        #Skew-T x-axis
        B = (x_in + (-y/p[0]))*p[0]
        x_out = (y_out+B)/p[0]

        return x_out, y_out

    #Process winds for the hodograph

def build_batch_lstsq_estimates(self, asset_returns, benchmark_returns):
        if not len(asset_returns) == len(benchmark_returns):
            raise '*WTF*'

        # Run Kalman filter on returns data
        beta = np.zeros(len(asset_returns))
        alpha = np.zeros(len(asset_returns))
        for enum_i, elem in enumerate(asset_returns):
            lookback = min(self.lookback, enum_i)
            # print '==> ', enum_i, len(asset_returns), len(beta)
            beta[enum_i], alpha[enum_i] = np.polyfit(benchmark_returns[enum_i - lookback:enum_i + 1],
                                                     asset_returns[enum_i - lookback:enum_i + 1], 1)

        # don't wanna do a line fit for less than 3 points, really
        beta[0], alpha[0] = 0, 0
        beta[1], alpha[1] = 0, 0

        self.alpha_betas_lstsq = np.array(zip(alpha, beta))

def early_stop_by_valid(valid_auc, k):
    '''Early stop"
    '''

    if k > len(valid_auc):
        return False
    print k, valid_auc

    x = range(k)
    y = valid_auc[-k:]

    # print x, y
    slop, bias = np.polyfit(x, y, 1)
    print "slop: {0:f}".format(slop)

    stopflag = False
    if slop <= FLAGS.slope:
        stopflag = True
    return stopflag

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def plot_scatter_charts(data, file_name):
    scatters = []
    for lang, values in data.items():
        s = figure(width=300, plot_height=300, title=lang)
        s.yaxis.formatter = NumeralTickFormatter(format="0.0a")
        s.circle(values[0], values[1], size=10, color="navy", alpha=0.5)
        x = np.linspace(1, 100, 10)
        # noinspection PyTupleAssignmentBalance
        m, b = np.polyfit(values[0], values[1], 1)
        y = m * x + b
        corr_coef = round(pearsonr(values[0], values[1])[0], 1)
        s.line(x, y, legend=f'PCC = {corr_coef}')
        scatters.append(s)
    split_scatters = split(scatters, 3)
    p = gridplot(split_scatters)
    output_file(file_name)
    show(p)

def defineRegionOfInterest(img, left_bottom = [0, 539], right_bottom = [900, 539], apex = [475, 320]):
    xsize = img.shape[1]
    ysize = img.shape[0]

    # Perform a linear fit (y=Ax+B) to each of the three sides of the triangle
    # np.polyfit returns the coefficients [A, B] of the fit
    fit_left = np.polyfit((left_bottom[0], apex[0]), (left_bottom[1], apex[1]), 1)
    fit_right = np.polyfit((right_bottom[0], apex[0]), (right_bottom[1], apex[1]), 1)
    fit_bottom = np.polyfit((left_bottom[0], right_bottom[0]), (left_bottom[1], right_bottom[1]), 1)

    # Find the region inside the lines
    XX, YY = np.meshgrid(np.arange(0, xsize), np.arange(0, ysize))

    region_thresholds = (YY > (XX*fit_left[0] + fit_left[1])) & \
                        (YY > (XX*fit_right[0] + fit_right[1])) & \
                        (YY < (XX*fit_bottom[0] + fit_bottom[1]))

    return region_thresholds

def testEncodeDecodeShift(self):
        x = np.linspace(-1, 1, 1000).astype(np.float32)
        with self.test_session() as sess:
            encoded = mu_law_encode(x, QUANT_LEVELS)
            decoded = mu_law_decode(encoded, QUANT_LEVELS)
            roundtripped = sess.run(decoded)

        # Detect non-unity scaling and non-zero shift in the roundtripped
        # signal by asserting that slope = 1 and y-intercept = 0 of line fit to
        # roundtripped vs x values.
        coeffs = np.polyfit(x, roundtripped, 1)
        slope = coeffs[0]
        y_intercept = coeffs[1]
        EPSILON = 1e-4
        self.assertNear(slope, 1.0, EPSILON)
        self.assertNear(y_intercept, 0.0, EPSILON)

def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

def hurst(ts):
    # Create the range of lag values
    window_len = 30
    hurst_ts = np.zeros((0,), dtype=np.int)
    for tail in range(window_len, len(ts)):
        lags = range(1, window_len // 2)
        # print(lags)
        # Calculate the array of the variances of the lagged differences
        cur_ts = ts[max(1, tail - window_len):tail]
        # cur_ts = ts[:tail]
        tau = [sqrt(std(subtract(cur_ts[lag:], cur_ts[:-lag]))) for lag in lags]

        # print("Time series slice len: ",len(cur_ts),len(tau),lags)
        # Use a linear fit to estimate the Hurst Exponent
        poly = polyfit(log(lags), log(tau), 1)

        # Return the Hurst exponent from the polyfit output
        # print(poly[0]*2.0)
        hurst_ts = np.append(hurst_ts, poly[0] * 2.0)
    return hurst_ts

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

def horner(circulation_time, times, temp):
    '''
    horner_bht_temp (circulation_time, times, temp)
    *Input parameters:
    - circulation_time - hours from last circulation;
    - times - total time since circulation stopped at 1st Run, 2nd Run and so on ...
    - temp - a list o temperatures coresponding to 1st Run, 2nd Run and so on ...
    *Returns:
    - horner_temp - formation temperature estimated by Horner method (thermometer readings
    from different runs)
    *Exemple of usage:
    horner(6, (7.0,11.5,19.5), (100,105,108))
        where:
        circulation_time = 6           # time since circulation stopped (hours)
        times = (7.0,11.5,19.5)        # total time since circulation stopped at 1st, 2nd, 3rd RUN (hours)
        temp=(100,105,108)             # temperatures recorded at 1st, 2nd, 3rd RUN (Celsius degrees)
    '''
    horner_time = np.array(times) / (circulation_time + np.array(times))
    slope,intercept = np.polyfit (np.log(horner_time), temp, 1)
    horner_temp=round(slope*np.log(1) +intercept,2)
    return horner_temp

def vander(points, deg):
    """N-dim Vandermonde matrix for data `points` and a polynomial of degree
    `deg`.

    Parameters
    ----------
    points : see polyfit()
    deg : int
        Degree of the poly (e.g. 3 for cubic).

    Returns
    -------
    vander : 2d array (npoints, (deg+1)**ndim)
    """
    powers = poly_powers(points.shape[1], deg)
    # low memory version, slower
    ##npoints = points.shape[0]
    ##vand = np.empty((npoints, (deg+1)**ndim), dtype=float)
    ##for ipoint in range(npoints):
    ##    vand[ipoint,:] = (points[ipoint]**powers).prod(axis=1)
    tmp = (points[...,None] ** np.swapaxes(powers, 0, 1)[None,...])
    return tmp.prod(axis=1)

def build_scatter_plot(dframe):
    #build scatter plot: stock1 vs. stock2
    daily_ret = get_daily_returns(dframe)
    symbols = []
    for symbol in dframe[0:]:
        symbols.append(symbol)
    daily_ret.plot(kind='scatter', x=symbols[0], y=symbols[1]) 
    #use ployfit() which needs x-coords and y-coords to fit a line -- y-coords are daily return values 
    #polyfit() returns first the polynomial coefficient (beta) and second the y intercept (alpha) -- in the form y = mx + b -- m is coefficient and b is intercept
    #the idea for plotting is for every value of x we find a value of y using the line equation y = mx +b
    beta_stock1, alpha_stock1 = np.polyfit(daily_ret[symbols[1]], daily_ret[symbols[0]], 1)
    beta_stock2, alpha_stock2 = np.polyfit(daily_ret[symbols[0]], daily_ret[symbols[1]], 1)

    print "beta of ", symbols[0], "= ", beta_stock1
    print "alpha of ", symbols[0], "= ", alpha_stock1   
    print "beta of ", symbols[1], "= ", beta_stock2
    print "alpha of ", symbols[1], "= ", alpha_stock2

    plt.plot(daily_ret[symbols[0]], beta_stock2 * daily_ret[symbols[0]] + alpha_stock2, '-', color='r')
    plt.show()

def fourierExtrapolation(x, n_predict):
    n = len(x)
    n_harm = 10                     # number of harmonics in model
    t = np.arange(0, n)
    p = np.polyfit(t, x, 1)         # find linear trend in x
    x_notrend = x - p[0] * t        # detrended x
    x_freqdom = fft.fft(x_notrend)  # detrended x in frequency domain
    f = fft.fftfreq(n)              # frequencies
    indexes = list(range(n))
    # sort indexes by frequency, lower -> higher
    indexes.sort(key = lambda i: np.absolute(f[i]))

    t = np.arange(0, n + n_predict)
    restored_sig = np.zeros(t.size)
    for i in indexes[:1 + n_harm * 2]:
        ampli = np.absolute(x_freqdom[i]) / n   # amplitude
        phase = np.angle(x_freqdom[i])          # phase
        restored_sig += ampli * np.cos(2 * np.pi * f[i] * t + phase)
    return restored_sig + p[0] * t