def test_polyfit(self):
     "Tests polyfit"
     # On ndarrays
     x = np.random.rand(10)
     y = np.random.rand(20).reshape(-1, 2)
     assert_almost_equal(polyfit(x, y, 3), np.polyfit(x, y, 3))
     # ON 1D maskedarrays
     x = x.view(MaskedArray)
     x[0] = masked
     y = y.view(MaskedArray)
     y[0, 0] = y[-1, -1] = masked
     #
     (C, R, K, S, D) = polyfit(x, y[:, 0], 3, full=True)
     (c, r, k, s, d) = np.polyfit(x[1:], y[1:, 0].compressed(), 3, full=True)
     for (a, a_) in zip((C, R, K, S, D), (c, r, k, s, d)):
         assert_almost_equal(a, a_)
     #
     (C, R, K, S, D) = polyfit(x, y[:, -1], 3, full=True)
     (c, r, k, s, d) = np.polyfit(x[1:-1], y[1:-1, -1], 3, full=True)
     for (a, a_) in zip((C, R, K, S, D), (c, r, k, s, d)):
         assert_almost_equal(a, a_)
     #
     (C, R, K, S, D) = polyfit(x, y, 3, full=True)
     (c, r, k, s, d) = np.polyfit(x[1:-1], y[1:-1, :], 3, full=True)
     for (a, a_) in zip((C, R, K, S, D), (c, r, k, s, d)):
         assert_almost_equal(a, a_)

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 meanclip3(xx,yy,slope, clipsig=3.0, maxiter=5, converge_num=0.1, verbose=0):
    from numpy import array, polyfit
    import numpy
    xx=array(xx)
    yy=array(yy)
    xx0=array(xx[:])
    yy0=array(yy[:])
    ct=len(yy)
    iter = 0; c1 = 1.0 ; c2 = 0.0
    while (c1 >= c2) and (iter < maxiter):
        lastct = ct
        pol = polyfit(xx0,yy0,1,full=True) ###
        mean0=pol[0][1]
        slope=pol[0][0]
        sig=numpy.std(yy0-mean0-slope*xx0)
        wsm = numpy.where( abs(yy0-xx0*slope) < mean0+clipsig*sig )
        ct = len(wsm[0])
        if ct > 0:
            xx0=xx0[wsm]
            yy0=yy0[wsm]
        c1 = abs(ct - lastct)
        c2 = converge_num * lastct
        iter += 1
# End of while loop
    pol = polyfit(xx0,yy0,1,full=True) ###
    mean0=pol[0][1]
    slope=pol[0][0]
    sig=numpy.std(yy0-mean0-slope*xx0)
    if verbose: pass
    return mean0, sig,slope,yy0,xx0

def test_run():
    # Read data
    dates = pd.date_range('2009-01-01', '2012-12-31')  # one month only
    symbols = ['SPY','XOM','GLD']
    df = get_data(symbols, dates)
    plot_data(df)

    # Compute daily returns
    daily_returns = compute_daily_returns(df)
    #plot_data(daily_returns, title="Daily returns", ylabel="Daily returns")

    # Scatterplot SPY vs XOM
    daily_returns.plot(kind='scatter',x='SPY',y='XOM')
    beta_XOM,alpha_XOM=np.polyfit(daily_returns['SPY'],daily_returns['XOM'],1)
    print "beta_XOM= ",beta_XOM
    print "alpha_XOM= ",alpha_XOM
    plt.plot(daily_returns['SPY'],beta_XOM*daily_returns['SPY']+alpha_XOM,'-',color='r')
    plt.grid()
    plt.show()

    # Scatterplot SPY vs GLD
    daily_returns.plot(kind='scatter',x='SPY',y='GLD')
    beta_GLD,alpha_GLD=np.polyfit(daily_returns['SPY'],daily_returns['GLD'],1)
    print "beta_GLD= ",beta_GLD
    print "alpha_GLD= ",alpha_GLD
    plt.plot(daily_returns['SPY'],beta_GLD*daily_returns['SPY']+alpha_GLD,'-',color='r')
    plt.grid()
    plt.show()

    # Calculate correlation coefficient
    print daily_returns.corr(method='pearson')

def polysmooth(x,y,z,NI,NJ):

    # size of the incoming array
    Nx, Ny = np.shape(z)
    x1d = x[:,0]
    y1d = y[0,:]

    # Get the C coefficients
    #NI = 7
    CIj = np.zeros((NI,Ny))
    for j in range (Ny):
        CIj[:,j] = np.flipud(np.polyfit(x1d,z[:,j],NI-1))

    # Get the D coefficients
    #NJ = 7
    DIJ = np.zeros((NI,NJ))
    for I in range (NI):
        DIJ[I,:] = np.flipud(np.polyfit(y1d,CIj[I,:],NJ-1))
    
    # Reconstruct the entire surface
    zsmooth = np.zeros((Nx,Ny))
    for I in range(NI):
        for J in range(NJ):
            zsmooth += DIJ[I,J]*x**I*y**J

    return zsmooth

def dfa(X, Ave = None, L = None, sampling= 1):
    """
    WIP on this function. It is basically copied and pasted from [PYEEG]_, without verification of the maths or unittests.
    """
    X = np.array(X)
    if Ave is None:
        Ave = np.mean(X)
    Y = np.cumsum(X)
    Y -= Ave
    if not L:
        max_power = np.int(np.log2(len(X)))-4
        L = X.size / 2 ** np.arange(4,max_power)
    if len(L)<2:
        raise Exception("Too few values for L. Time series too short?")
    F = np.zeros(len(L)) # F(n) of different given box length n

    for i,n in enumerate(L):
        sampled = 0
        for j in xrange(0,len(X) -n ,n):

            if np.random.rand() < sampling:
                F[i] += np.polyfit(np.arange(j,j+n), Y[j:j+n],1, full=True)[1]
                sampled += 1
        if sampled > 0:
            F[i] /= float(sampled)

    LF = np.array([(l,f) for l,f in zip(L,F) if l>0]).T

    F = np.sqrt(LF[1])
    Alpha = np.polyfit(np.log(LF[0]), np.log(F),1)[0]
    return Alpha

def fit(data, nz):
    x = [0 for iz in range(0, nz, 1)]
    y = [0 for iz in range(0, nz, 1)]
    z = [iz for iz in range(0, nz, 1)]

    for iz in range(0, nz, 1):
        x[iz], y[iz] = ndimage.measurements.center_of_mass(np.array(data[:,:,iz]))

    #Fit centerline in the Z-X plane using polynomial function
    print '\nFit centerline in the Z-X plane using polynomial function...'
    coeffsx = np.polyfit(z, x, 1)
    polyx = np.poly1d(coeffsx)
    x_fit = np.polyval(polyx, z)
    print 'x_fit'
    print x_fit

    #Fit centerline in the Z-Y plane using polynomial function
    print '\nFit centerline in the Z-Y plane using polynomial function...'
    coeffsy = np.polyfit(z, y, 1)
    polyy = np.poly1d(coeffsy)
    y_fit = np.polyval(polyy, z)


    #### 3D plot
    fig1 = plt.figure()
    ax = Axes3D(fig1)
    ax.plot(x,y,z,zdir='z')
    ax.plot(x_fit,y_fit,z,zdir='z')
    plt.show()
    return x, y, x_fit, y_fit

def fitPlotFW(nGraphs, nNodesIni, nNodesFin, step,sparseFactor =0.25):

    times1 = timeDijkstraMAllPairs(nGraphs, nNodesIni, nNodesFin, step, sparseFactor)
    times2 = timeFloydWarshall(nGraphs, nNodesIni, nNodesFin, step, sparseFactor)

    x = []
    for z in range(nNodesIni, nNodesFin, step):
        x.append(z)

    ydata1 = np.polyfit(x, times1, 2)
    f1 = np.poly1d(ydata1)
    yfit1 = np.linspace(nNodesIni, nNodesFin,nNodesFin - nNodesIni / step)

    ydata2 = np.polyfit(x, times2, 2)
    f2 = np.poly1d(ydata2)
    yfit2 = np.linspace(nNodesIni, nNodesFin,nNodesFin - nNodesIni / step)

    plt.plot(x,times1,"b.", label="valor real Dijkstra")
    plt.plot(yfit1, f1(yfit1),"r-", label='Dijkstra')
    plt.plot(x,times2, "r.", label="valor real FloydWarshall")
    plt.plot(yfit2,f2(yfit2),"b-", label ='FloydWarshall')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    return

def poleward_speed(KD):

    N = int(KD.params['N'])
    trans_MA = KD.params['t_A']
    dt = KD.params['dt']
    start = int(trans_MA/dt)

    trans_MA = KD.params['t_A']
    dt = KD.params['dt']
    trans_AB = anaph_transAB(KD)

    pole_speeds = []
    for ch in KD.chromosomes:
        r_stop = int(ch.cen_A.toa/dt)
        if r_stop - start > 2:
            r_dist = - KD.spbR.traj[start:r_stop] + ch.cen_A.traj[start:r_stop]
            elapsed = np.r_[trans_MA:ch.cen_A.toa:dt]
            (ra,rb) = np.polyfit(elapsed, r_dist, 1)
            pole_speeds.append(ra)
        l_stop = int(ch.cen_B.toa/dt)
        if l_stop - start > 2:
            l_dist = KD.spbL.traj[start:l_stop] - ch.cen_B.traj[start:l_stop]
            elapsed = np.r_[trans_MA:ch.cen_B.toa:dt]
            (la,lb) = np.polyfit(elapsed, l_dist, 1)
            pole_speeds.append(la)

    pole_speeds = np.array(pole_speeds)

    return pole_speeds.mean(), pole_speeds.std()

    def FitLine(self, definition):
        '''
        Fits a 1st order polynom (line) to the
        start and end points of a straight path,
        then populates the line with equally spaced
        points, and returns the list of points
        '''
    
        Ax = self.A.get_Pos_X();
        Ay = self.A.get_Pos_Y();
        
        Bx = self.B.get_Pos_X();
        By = self.B.get_Pos_Y();
        
        SubWP_No = numpy.linalg.norm(numpy.array([Ax-Bx,Ay-By])) * definition * 0.01;
        '''
        If the path is vertical, the X and Y axes must be swapped before the
        polynom fitting and populating, then switched back to return the
        proper point coordinates
        ''' 
        if abs(Ax - Bx) < 1:

            self.poly = numpy.polyfit([Ay, By], [Ax, Bx], 1);
            prange = numpy.linspace(Ay, By, SubWP_No);

            values = numpy.polyval(self.poly, prange);
            self.SubWP = numpy.array([prange, values]);
            
        else:            
            self.poly = numpy.polyfit([Ax, Bx], [Ay, By], 1);
            prange = numpy.linspace(Ax, Bx, SubWP_No);
            values = numpy.polyval(self.poly, prange);
            self.SubWP = numpy.array([values, prange]);

def polyfit(x, y, deg, rcond=None, full=False):
    'units wrapped numpy.polyfit'
    
    _polyfit = np.polyfit

    if full:
        dP, residuals, rank, singular_values, rcond  = np.polyfit(np.array(x), np.array(y), deg, rcond, full)
    else:
        dP = np.polyfit(np.array(x), np.array(y), deg, rcond, full)

    # now we need to put units on P
    # p(x) = p[0] * x**deg + ... + p[deg]
    P = []
    for i, p in enumerate(dP):
        power = deg - i
        X = x**power

        # units on ith element of P from y / X
        uX = Unit(1.0, X.exponents, X.label)
        uy = Unit(1.0, y.exponents, y.label)
        uPi = uy / uX
        # so annoying. if you do dP[i] * uP you lose units.
        P += [uPi * p]

    if full:
        return P, residuals, rank, singular_values, rcond
    else:
        return P

 def homogenize(self, other, maxdiff=1):
     """ Return overlapping part of self and other as (self, other) tuple.
     Homogenize intensities so that the images can be used with
     combine_frequencies. Note that this works best when most of the 
     picture is signal, so use :py:meth:`in_interval` to select the subset
     of your image before applying this method.
     
     Parameters
     ----------
     other : CallistoSpectrogram
         Spectrogram to be homogenized with the current one.
     maxdiff : float
         Threshold for which frequencies are considered equal.
     """
     one, two = self._overlap(other)
     pairs_indices, factors, constants = one._homogenize_params(
         two, maxdiff
     )
     # XXX: Maybe (xd.freq_axis[x] + yd.freq_axis[y]) / 2.
     pairs_freqs = [one.freq_axis[x] for x, y in pairs_indices]
     
     # XXX: Extrapolation does not work this way.
     # XXX: Improve.
     f1 = np.polyfit(pairs_freqs, factors, 3)
     f2 = np.polyfit(pairs_freqs, constants, 3)
     
     return (
         one,
         two * polyfun_at(f1, two.freq_axis)[:, np.newaxis] +
             polyfun_at(f2, two.freq_axis)[:, np.newaxis]
     )

def plot_reml(C_eta, sel, sel_label, ax, ax2, sigma_diag):
    coh_arr = np.linspace(0.00, 0.30, 16)
    reml_arr = reml(C_eta, coh_arr, sel, sel_label, ax2, sigma_diag)
    ax.plot(coh_arr, reml_arr, 'm+', label=sel_label)
    ax.set_xlabel(r'$\sigma_i$')
    ax.set_ylabel('REML')
    ax.set_xlim(-0.02,0.30)

    # Fit a polynomial to the points
    coeff=np.polyfit(coh_arr, reml_arr, 4)
    p=np.poly1d(coeff)
    coh_arr2=np.linspace(0.00, 0.20, 201)
    fit=p(coh_arr2)
    ax.plot(coh_arr2,fit)

    # Determine where the minumum occurs
    minimum=sp.fmin(p,0.1,disp=False)
    print "Miminum at %4.3f" % (minimum[0])
    # The uncetainty occurs whwn the the REML increases by 1 - need to double check
    # To find where the the REML increases by 1, we look for the roots
    coeff=np.polyfit(coh_arr, reml_arr-p(minimum[0])-1, 4)
    p=np.poly1d(coeff)
    sol=sp.root(p,[0.0,0.2])

    m=minimum[0]
    upper=sol.x[1]-m
    lower=m-sol.x[0]
    ax.plot(coh_arr2,fit,label="Min at %4.2f+%4.2f-%4.2f" % (m,upper,lower))
    
    return

    def fit(self, data):
        magnitude = data[0]
        time = data[1]

        global m_21
        global m_31
        global m_32

        Nsf = 100
        Np = 100
        sf1 = np.zeros(Nsf)
        sf2 = np.zeros(Nsf)
        sf3 = np.zeros(Nsf)
        f = interp1d(time, magnitude)

        time_int = np.linspace(np.min(time), np.max(time), Np)
        mag_int = f(time_int)

        for tau in np.arange(1, Nsf):
            sf1[tau-1] = np.mean(np.power(np.abs(mag_int[0:Np-tau] - mag_int[tau:Np]) , 1.0))
            sf2[tau-1] = np.mean(np.abs(np.power(np.abs(mag_int[0:Np-tau] - mag_int[tau:Np]) , 2.0)))
            sf3[tau-1] = np.mean(np.abs(np.power(np.abs(mag_int[0:Np-tau] - mag_int[tau:Np]) , 3.0)))
        sf1_log = np.log10(np.trim_zeros(sf1))
        sf2_log = np.log10(np.trim_zeros(sf2))
        sf3_log = np.log10(np.trim_zeros(sf3))

        m_21, b_21 = np.polyfit(sf1_log, sf2_log, 1)
        m_31, b_31 = np.polyfit(sf1_log, sf3_log, 1)
        m_32, b_32 = np.polyfit(sf2_log, sf3_log, 1)

        return m_21

def calc_def_potential(info_file):
	data =[]
	counter = 0
	with open(info_file, 'r') as energy_info:
#		enfo.seek(0)
#		next(enfo)
			for row in energy_info:
				counter += 1
				if counter > 1:
					data.append([float(i) for i in row.split()])
	
	data = np.mat(data)
	eqidx = np.where(data[:,0]==0)[0]
#	eqidx = abs(int(data[0,0])) + 1
	eqvol = data[eqidx,1]
	data[:,1] -= eqvol 
	data[:,4] -= data[:,3]
	data[:,5] -= data[:,3]
	
	x = np.array(np.ravel(data[:,1]))
	y = np.array(np.ravel(data[:,4]))
	z = np.polyfit(x, y, 6)
	VBM_deformation = abs(z[-2]*eqvol)
	
	x = np.array(np.ravel(data[:,1]))
	y = np.array(np.ravel(data[:,5]))
	z = np.polyfit(x, y, 6)
	CBM_deformation = abs(z[-2]*eqvol)
	
	return float(VBM_deformation), float(CBM_deformation)

    def get_spectrum(self,sample_str,x_scale = 'energy'):

        #normalize the signals to the tube current
        direct_beam = self.scan_groups['direct_beam']
        sample = self.scan_groups[sample_str]

        theta = direct_beam['signal']['theta']

        I0 = direct_beam['signal']['intensity']/direct_beam['signal']['tube_current']
        I0_err = direct_beam['signal']['intensity_error']/direct_beam['signal']['tube_current']

        I = sample['signal']['intensity']/sample['signal']['tube_current']
        I_err = sample['signal']['intensity_error']/sample['signal']['tube_current']

        theta_I0bg = direct_beam['background']['theta']
        theta_Ibg = sample['background']['theta']

        I0_bg =  direct_beam['background']['intensity']/direct_beam['background']['tube_current']
        I_bg =  sample['background']['intensity']/sample['background']['tube_current']

        #fit backgrounds
        p0 = np.polyfit(theta_I0bg,I0_bg,self.background_fit_order)
        p = np.polyfit(theta_Ibg,I_bg,self.background_fit_order)

        #compute mux
        mux = -np.log((I-np.polyval(p,theta))/(I0-np.polyval(p0,theta)))
        mux_error = np.sqrt((I0_err/I0)**2 + (I_err/I)**2)

        if x_scale == 'theta':
            return theta+self.theta_calibration, mux, mux_error
        else:
            return energy(theta+self.theta_calibration,*self.analyser), mux, mux_error

def DFA(indata,scale,q,m):
    y = np.cumsum(indata-indata.mean())             #Equation 1 in paper
    RMSt = []                                       #Temporary RMS variable: contain F(s,v) value
    F = []                                          #F: Fluctuation function
    N = len(indata)
    print 'len indata: ',N
    for i in range(len(scale)):
        ns = int(np.floor(len(y)/scale[i]))         #number of segments: Ns = int(N/s)
        for v in range(2*ns):
            if v < ns:
                index_start = v*scale[i]
                index_end = (v+1)*scale[i]
            else:
                index_start = N - (v-ns)*scale[i]-scale[i]
                index_end = N - (v-ns)*scale[i]
            index = range(index_start,index_end)    #calculate index for each segment
            yv = y[index_start:index_end]           #Extract values of time series for each segments
            c = np.polyfit(index,yv,m)
            fit = np.polyval(c,index)
            RMSt.append(math.sqrt(np.mean((yv-fit)**2))) #Equation 2. But calculating only F(v,s) not F(v,s)**2
        RMS = np.asarray(RMSt)                      #Convert RMSt to array
        qRMS = RMS**q
        F.append(np.mean(qRMS)**(1.0/q))              #Equation 4
        del RMSt[:]                                 #Reset RMSt[:]
    C = np.polyfit(np.log2(scale),np.log2(F),1)
    H = C[0]                                        #Hurst parameter
    return (H,scale,F)

def optimum_polyfit(x, y, score=functoolz.compose(np.max, np.abs), max_degree=50, stop_at=1e-10):
    """
    Optimize the degree of a polyfit polynomial so that score(y - poly(x)) is minimized.

    :param max_degree: The maximum degree to try. LinAlgErrors are automatically ignored.
    :param stop_at: If a score lower than this is reached, the function returns early
    :param score: The score function that is applied to y - poly(x). Default: max deviation.
    :return A tuple (poly1d object, degree, score)
    """
    scores = np.empty(max_degree - 1, dtype=np.float64)
    # Ignore rank warnings now, but do not ignore for the final polynomial if not early returning
    with warnings.catch_warnings():
        warnings.simplefilter('ignore', np.RankWarning)
        for deg in range(1, max_degree):
            # Set score to max float value
            try:
                poly = np.poly1d(np.polyfit(x, y, deg))
            except np.linalg.LinAlgError:
                scores[deg - 1] = np.finfo(np.float64).max
                continue
            scores[deg - 1] = score(y - poly(x))
            # Early return if we found a polynomial that is good enough
            if scores[deg - 1] <= stop_at:
                return poly, deg, scores[deg - 1]
    # Find minimum score
    deg = np.argmin(scores) + 1
    # Compute polyfit for that degreet
    poly = np.poly1d(np.polyfit(x, y, deg))
    return poly, deg, np.min(scores)

def next_start(X, F, G, i, j,scandirec=0,radius=1.,method='default'):
        radius*=np.sqrt(abs(F[i-1]**2+G[i-1]**2))
        df=j*radius*np.cos(scandirec)
        dg=j*radius*np.sin(scandirec)
        if 'default' in method: #use old root position as starting point, plus potential scandirec modifications
                fval=F[i]+df
                gval=G[i]+dg
        elif 'predict' in method: #use first order prediction of the root's new position
                if i>1:
                        p=np.polyfit(X[i-2:i+1],F[i-2:i+1],2)
                        fval=p[0]*X[i+1]**2+p[1]*X[i+1]+p[2]+df
                        p=np.polyfit(X[i-2:i+1],G[i-2:i+1],2)
                        gval=p[0]*X[i+1]**2+p[1]*X[i+1]+p[2]+dg                     
                elif i==1:
                        fval=F[i]+np.diff(F)[i-1]/np.diff(X)[i-1]*np.diff(X)[i]+df
                        gval=G[i]+np.diff(G)[i-1]/np.diff(X)[i-1]*np.diff(X)[i]+dg
                else:
                        fval=F[i]
                        gval=G[i]
                if 'ldr' in method:
                        gval=-0.75*abs(G[i]+np.diff(G)[i-1]/np.diff(X)[i-1]*np.diff(X)[i])+dg
        elif 'ldr' in method: #lower damping rate (designed to prefer more weakly damped modes)
                gval=-0.75*abs(G[i])+dg #always aim for lower damping rate so the less damped mode will be chosen
                fval=F[i]+0.03*F[i]*np.sign(np.diff(X)[-1])+df #include slight frequency shift (depending on the scan direction) in order to help jump across resonances
        #lower limit for damping rates, to prevent runaway to extremely small values 
        if abs(gval)<1e-13: gval=np.sign(gval)*1e-13
        return [fval,gval] 

    def estimate_linear_fit(self, data, split_b, less_than=True):
        """Estimate a linear fit by taking log of data.

        Parameters
        ----------
        data : array
            An array containing the data to be fit

        split_b : float
            The b value to split the data

        less_than : bool
            If True, splitting occurs for bvalues less than split_b

        Returns
        -------
        S0 : float
            The estimated S0 value. (intercept)

        D : float
            The estimated value of D.
        """
        if less_than:
            bvals_split = self.gtab.bvals[self.gtab.bvals <= split_b]
            D, neg_log_S0 = np.polyfit(bvals_split,
                                       -np.log(data[self.gtab.bvals <=
                                                    split_b]), 1)
        else:
            bvals_split = self.gtab.bvals[self.gtab.bvals >= split_b]
            D, neg_log_S0 = np.polyfit(bvals_split,
                                       -np.log(data[self.gtab.bvals >=
                                                    split_b]), 1)

        S0 = np.exp(-neg_log_S0)
        return S0, D