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

def polynom_approximate(x, y, polynom_degree):
    # Get polynom params
    polynom = sp.polyfit(x, y, polynom_degree)

    return polynom

# Model error calculate

def fit_linear_model(x,y,slope=None):
    '''Linear least squares (LSQ) linear regression (polynomial fit of degree=1)
    Returns:
        m (float) = slope of linear regression line of form (y = m*x + b)
        b (float) = intercept of linear regression line'''

    assert len(x)==len(y), ("Arrays x & Y must be equal length to fit "
                            "linear regression model.")
    if slope == None:
        (m,b) = scipy.polyfit(x,y,deg=1)
    else:
        LSQ = lambda b: np.sum( (y-(slope*x+b))**2.0 )
        res = scipy.optimize.minimize(LSQ,x0=1,bounds=None)
        (m,b) = (slope,res.x[0])
    return (m,b)

def bowtie_polynom(modis_img,cs,folder):
    print 'Determine overlap pattern... '
    sw=10000/cs #stripwidth
    overlaplist=[]#define list to store number of overlapped lines
    #devide in parts with a width of 40 pixel
    for i in sp.arange(0,modis_img.shape[1]-40,40):
        part=modis_img[:,i:i+39]
        #search in every scanning strip
        samples=[]
        for j in sp.arange(sw-2,part.shape[0]-sw,sw):
            target=part[j-1:j+1,:] #cut out a target, which overlapped counter-part shall be found
            searchwindow=part[j+2:j+sw+2] #,: cut out the window, where the overlapped counter part might be located
            #start the search
            c=[] #calculate correlation coefficients of every given offset from 3 to 11
            for offset in sp.arange(3,sw/2+1):
                imgpart=searchwindow[offset-3:offset-1] #,: cut out image, which has to be compared with the target
                c.append(sp.corrcoef(imgpart.flatten(),target.flatten())[0,1])#calculate correlatoin coefficient
            c=sp.array(c)
            overl=sp.ndimage.measurements.maximum_position(c)[0]+3 #find the overlap with the highes correlation coefficient
            samples.append([overl,c.max()]) #attach overlap and correlation coefficient to the sample list
        samples=sp.array(samples)
        #print i, samples[:,1].mean()
        if samples[:,1].mean() > 0.9: #chek the mean correlation coefficient:
            #print('Bowtie Correlation high - removing effect')
            overlaplist.append([i+20,samples[:,0].mean()]) #save result, if correlation coefficient is high
            #print(overlaplist)
            o=sp.array(overlaplist)
            X=o[:,0]
            overlap=o[:,1]
            #Calculate a second order Polynom to describe the overlap
            p=sp.polyfit(X,overlap,2)
            #print 'done, Overlap polynom: '+str(p)
        else:
            #print('low Bowtie correlation')
            p = [1.,  1.,  1.]
            #overlaplist.append([i+20,1])
            #os.system('rm -r '+folder)
            #print('scene deleted') 
    return p

def interp_helper(all_data, trend_data, time_from):
    'performs lf spline + hf fft interpolation of radial distance'
    all_times, all_values = zip(*all_data)
    trend_times, trend_values = zip(*trend_data)

    split_time = int(time_to_index(time_from, all_times[0]))

    trend_indices = array([time_to_index(item, all_times[0]) for item in trend_times])
    spline = splrep(trend_indices, array(trend_values))

    all_indices = array([time_to_index(item, all_times[0]) for item in all_times])
    trend = splev(all_indices, spline)
    detrended = array(all_values) - trend
    trend_add = splev(arange(split_time, all_indices[-1]+1), spline)

    dense_samples = detrended[:split_time]
    sparse_samples = detrended[split_time:]
    sparse_indices = (all_indices[split_time:]-split_time).astype(int)
    amp = log(absolute(rfft(dense_samples)))
    dense_freq = rfftfreq(dense_samples.size, 5)
    periods = (3000.0, 300.0)
    ind_from = int(round(1/(periods[0]*dense_freq[1])))
    ind_to = int(round(1/(periods[1]*dense_freq[1])))
    slope, _ = polyfit(log(dense_freq[ind_from:ind_to]), amp[ind_from:ind_to], 1)

    params = {
        't_max': periods[0],
        'slope': slope,
        'n_harm': 9,
        'scale': [20, 4, 2*pi]
    }
    series_func, residual_func = make_residual_func(sparse_samples, sparse_indices, **params)

    x0 = array([0.5]*(params["n_harm"]+2))
    bounds = [(0, 1)]*(params["n_harm"]+2)
    result = minimize(residual_func, x0, method="L-BFGS-B", bounds=bounds, options={'eps':1e-2})
    interp_values = [trend + high_freq for trend, high_freq in
                     zip(trend_add, series_func(result.x)[:sparse_indices[-1]+1])]
    #make_qc_plot(arange(sparse_indices[-1]+1), interp_values,
    #             sparse_indices, array(all_values[split_time:]))
    interp_times = [index_to_time(ind, time_from) for ind in range(sparse_indices[-1]+1)]
    return list(zip(interp_times, interp_values))