def lagrange(self):
        '''Lagrange Interpolation. Unsuccessful because it's a linear equation, so doesn't
        get the target result that I want. The result is -3565+2421.41*(x)'''

        orig_stdout = sys.stdout
        f = file('lagrangeout.txt', 'w')
        sys.stdout = f

        tmp = scipy.poly1d([0])
        result=scipy.poly1d([0])
        for value,key in zip(self.complaint_allstate.values(),self.complaint_allstate.keys()):
            for i in value.keys():
                numerator=scipy.poly1d([1])
                denom = 1.0
                for j in value.keys():
                    if (i != j):
                        tmp = scipy.poly1d([1,-j])
                        numerator = numerator * tmp
                        denom = denom * (i - j)
                tmp = (numerator/denom) * value.get(i)
                result = result + tmp
            print key + " : "
            print result
            print '\n'
        sys.stdout = orig_stdout
        f.close()

    def trainingAndTesting(self):
        global train
        global x, y, xa, xb, ya, yb
        # separating training from testing data
        frac = 0.3
        split_idx = int(frac * len(xb))
        rangeX = range(len(xb))
        listX = list(rangeX)
        logging.info("delta : %i", len(set(rangeX).difference(listX)))

        shuffled = sp.random.permutation(list(range(len(xb))))
        test = sorted(shuffled[:split_idx])

        train = sorted(shuffled[split_idx:])
        fbt1 = sp.poly1d(sp.polyfit(xb[train], yb[train], 1))
        fbt2 = sp.poly1d(sp.polyfit(xb[train], yb[train], 2))
        print("fbt2(x)= \n%s" % fbt2)
        print("fbt2(x)-100,000= \n%s" % (fbt2 - 100000))
        print("fbt2(x)= \n%s" % (fbt2))
        fbt3 = sp.poly1d(sp.polyfit(xb[train], yb[train], 3))
        fbt10 = sp.poly1d(sp.polyfit(xb[train], yb[train], 10))
        fbt100 = sp.poly1d(sp.polyfit(xb[train], yb[train], 100))

        print("Test errors for only the time after inflection point")
        for f in [fbt1, fbt2, fbt3, fbt10, fbt100]:
            print("Error d=%i: %f" % (f.order, self.error(f, xb[test], yb[test])))

def main():
    data = sp.genfromtxt('./data/web_traffic.tsv', delimiter='\t')
    x = data[:, 0]
    y = data[:, 1]
    x = x[~sp.isnan(y)]
    y = y[~sp.isnan(y)]
    fp1 = sp.polyfit(x, y, 1)
    print('Model parameters for fp1 %s' % fp1)
    f1 = sp.poly1d(fp1)
    print('This is the error rate for fp1 %f' % error(f1, x, y))

    fp2 = sp.polyfit(x, y, 2)
    print('Model parameters for fp2 %s' % fp2)
    f2 = sp.poly1d(fp2)
    print('This is the error rate for fp2 %f' % error(f2, x, y))

    plt.scatter(x, y,color= 'pink')
    plt.title('My first impression')
    plt.xlabel('Time')
    plt.ylabel('#Hits')
    plt.xticks([w * 7 * 24 for w in range(10)], ['week %i' % w for w in range(10)])
    fx = sp.linspace(0, x[-1], 1000)
    plt.plot(fx, f1(fx), linewidth=3,color='cyan')


    plt.plot(fx, f2(fx), linewidth=3, linestyle='--',color= 'red')
    plt.legend(['d = %i' %f1.order, 'd = %i' %f2.order], loc='upper left')
    plt.autoscale(tight=True)
    plt.grid()
    plt.show()

def find_roots_simple(B=0.0, kr=0.0, kz=0.0, Omega=0.0, zeta=0.0, m=1.0,
                      nu=2.98e-3, eta=2.57e3, rho=6.36, r1=7.06, r2=20.3):
    """Solve the local dispersion relation."""
    pi = math.pi

    ktheta = m/(0.5*(r2 + r1))
    k = sqrt(kr**2 + kz**2 + ktheta**2)
    va = B*1.0/math.sqrt(4.0*pi*rho)
    wa = kz*va #Alfven frequency \omega_A

    #gn = \gamma_n = \gamma + \eta k^2 = i\omega + \eta k**2
    gn = scipy.poly1d([1, eta*k**2])
    gv = scipy.poly1d([1, nu*k**2])
    
    dr = (gn*gv + wa**2)**2 + 2*Omega**2*zeta*(kz**2/k**2)*gn**2 \
         + 2*Omega**2*(zeta-2)*(kz**2/k**2)*wa**2
    
    dr_roots = roots(dr)
    
    #Convert to \omega
    #Since \gamma t = -i \omega t, \omega = i \gamma
    
    dr_roots_return = sort(1j*dr_roots)
    
    return dr_roots_return

    def fitKvsPower(self, filename="PowerVGain.dat", zeroed="_0"):

        if not hasattr(self, "k_avg"):
            raise RuntimeError, "Must load all data and run calculate() before fitting"

        gain, self.power = loadtxt(filename, skiprows=1, unpack=True)

        if zeroed != None:
            gain_0, power_0 = loadtxt(filename.replace(".", zeroed + ".", 1), skiprows=1, unpack=True)
            self.power = self.power - power_0
            savetxt(filename.replace(".", "_subtract."), self.power, fmt="%.5f")

        if self.gain.tolist() != gain.tolist():
            raise ValueError, "Laser power was not measured over the same gains as calibration measurements:\n\
		Laser gain: " + str(
                gain
            ) + "\nMeasurement gains: " + str(
                self.gain
            )

        self.kfitX = polyfit(self.power, self.k_avg[:, 0], 1)
        self.kfitY = polyfit(self.power, self.k_avg[:, 1], 1)
        print "Fit parameters for X axis:"
        print poly1d(self.kfitX, variable="power")
        print "\nFit parameters for Y axis:"
        print poly1d(self.kfitY, variable="power")

def _systopoly1d(sys):
    """Extract numerator and denominator polynomails for a system"""
    # Allow inputs from the signal processing toolbox
    if (isinstance(sys, scipy.signal.lti)):
        nump = sys.num
        denp = sys.den

    else:
        # Convert to a transfer function, if needed
        sys = _convert_to_transfer_function(sys)

        # Make sure we have a SISO system
        if (sys.inputs > 1 or sys.outputs > 1):
            raise ControlMIMONotImplemented()

        # Start by extracting the numerator and denominator from system object
        nump = sys.num[0][0]
        denp = sys.den[0][0]

    # Check to see if num, den are already polynomials; otherwise convert
    if (not isinstance(nump, poly1d)):
        nump = poly1d(nump)

    if (not isinstance(denp, poly1d)):
        denp = poly1d(denp)

    return (nump, denp)

def plotFullModel(dm, sacc):
	
	"""
	"""


	fig = plt.figure(figsize = (15, 5))
	plt.suptitle("Sacc = %s" % sacc)
	ax1 = plt.subplot(141)
	lCols = ["blue", "red", "orange", "yellow", "green", "pink"]
	
	xVar = "sacc%s_ex" % sacc
	yVar = "sacc%s_ey" % sacc
	
	dm = dm.select("%s != ''" % xVar)
	dm = dm.select("%s != ''" % yVar)

	dm = dm.select("%s != -1000" % xVar)
	dm = dm.select("%s != -1000" % yVar)

	for a in dm.unique("realAngle"):
		_dm = dm.select("realAngle == %s" % a)
		col = lCols.pop()
		plt.scatter(_dm[xVar], _dm[yVar], color = col, marker = ".", label = a)
	plt.axvline(constants.xCen, color = gray[3], linestyle = '--')
	plt.axhline(constants.yCen, color = gray[3], linestyle = '--')

	ax2 = plt.subplot(142)
	pm = PivotMatrix(dm, ["stim_type", "realAngle"], ["file"], "xNorm1", colsWithin =True)
	pm.linePlot(fig = fig)
	pm.save("PM.csv")
	pm._print()
	plt.axhline(0, color = gray[3], linestyle = "--")

	ax3 = plt.subplot(143)
	plt.title("object")
	dmObj= dm.select("stim_type == 'object'")
	#slope, intercept, r_value, p_value, std_err  = scipy.stats.linregress(stimObj["ecc"], stimObj["xNorm1"])
	x = dmObj["ecc"]
	y = dmObj["xNorm%s" % sacc]
	fit = scipy.polyfit(x,y,1)
	fit_fn = scipy.poly1d(fit)
	plt.plot(x,y, 'yo', x, fit_fn(x), '--k')

	ax4 = plt.subplot(144)
	plt.title("non-object")
	dmNo= dm.select("stim_type == 'non-object'")
	#slope, intercept, r_value, p_value, std_err  = scipy.stats.linregress(stimObj["ecc"], stimObj["xNorm1"])
	x = dmNo["ecc"]
	y = dmNo["xNorm%s" % sacc]
	fit = scipy.polyfit(x,y,1)
	fit_fn = scipy.poly1d(fit)
	plt.plot(x,y, 'yo', x, fit_fn(x), '--k')

	plt.savefig("./plots/Full_model_004C_sacc%s.png" % sacc)
	plt.show()

def find_roots(B=0.0, kr=0.0, kz=0.0, Omega=0.0, zeta=0.0, m=0,
               nu=2.98e-3, eta=2.57e3, rho=6.36, r1=7.06, r2=20.3):

    #Parameters are in cgs units (Gauss, cm, 1/cm, cm^2/sec, gm/cm^3, etc)
    #We will be solving for the Doppler shifted frequency
    #\bar{\omega} = \omega - m\Omega
    
    pi = math.pi
    
    ktheta = m/(0.5*(r2+r1))
    k = sqrt(kr**2 + kz**2 + ktheta**2)
    va = B*1.0/math.sqrt(4.0*pi*rho)
    wa = kz*va #Alfven frequency \omega_A
    wr = (zeta-2)*Omega*kr*ktheta/k**2 #Rossby wave frequency \omega_R
    kfrac = kz/k
    
    #wn = \omega_\eta = \bar{\omega} - i\eta k^2
    #wv = \omega_\nu = \bar{\omega} - i\nu k^2
    
    wn = scipy.poly1d([1, -1j*eta*k**2])
    wv = scipy.poly1d([1, -1j*nu*k**2])

    #Note that for all of these terms are start off with the operations
    #on the poly1d objects, because otherwise I think I've seen instances
    #where the operator overloading fails, and multiplies, adds, and
    #exponentiations don't treat the objects correctly.

    #Term 1: \omega_{\eta}(\omega_A^2 - \omega_{\eta}\omega_{\nu})^2
    
    term1 = wn*(-wn*wv + wa**2)**2
    
    #Term 2: -2\zeta\Omega^2(k_z^2/k^2)\omega_{\eta}^3
    
    term2 = -wn**3*2.0*zeta*Omega**2*kfrac**2
    
    #Term 3: 2(\zeta-2)\Omega^2\omega_A^2(k_z^2/k^2)\omega_{\eta}
    
    term3 = wn*2.0*(zeta-2)*Omega**2*wa**2*kfrac**2
    
    #Term 4: i*\omega_R*(\omega_A^2 - \omega_{\nu}\omega_{\eta})*
    #(\omega_A^2 - \omega_{\eta}^2)
    #where \omega_R = (\zeta-2)\Omega k_r k_{\theta} / k^2
    
    term4 = (-wv*wn + wa**2)*(-wn**2 + wa**2)*1j*wr
    
    dr = term1 + term2 + term3 + term4

    #Minus sign here, because this whole thing was derived with
    #a bizarre assumed dependence as exp(\omega t - k\cdot x), when
    #a reasonable right-going wave is exp(k\cdot x - \omega t).
    
    dr_roots = -roots(dr)

    return sort(dr_roots)

def lagrangeInterpolation (x, gx):
   #setting result = 0
   result = scipy.poly1d([0.0]) 
   for i in range(0,len(x)): #number of polynomials L_k(x).
      temp_numerator = scipy.poly1d([1.0]) # resets temp_numerator such that a new numerator can be created for each i.
      denumerator = 1.0 #resets denumerator such that a new denumerator can be created for each i.
      for j in range(0,len(x)):
          if i != j:
              temp_numerator *= scipy.poly1d([1.0,-x[j]]) #finds numerator for L_i
              denumerator *= x[i]-x[j] #finds denumerator for L_i
      result += (temp_numerator/denumerator) * gx[i] #linear combination
   return result;

    def save(self, filename='calibration.dat', delimiter='\t', raw=True):

	# save calculated values in comments
	comments=\
	    """# ky vs power: %s # kx vs power: %s\n""" % \
		(str(poly1d(self.kfitX, variable='power')).strip(), 
		str(poly1d(self.kfitY, variable='power')).strip() )

	# patterns for col names of raw data
	lro_raw_headings = ['kxVar%d', 'kyVar%d', 'kxLRO%d', 'kyLRO%d']
	var_raw_headings = ['kx%d', 'ky%d']

	# col names for raw data
	if raw:
	    apply_to = lambda mylist, value: [ x % value for x in mylist ]
	    lro_head = []
	    var_head = []
	    for lronum, varnum in zip(self.__lro_file_list, self.__var_file_list):
		lro_head += apply_to(lro_raw_headings, lronum)
		var_head += apply_to(var_raw_headings, varnum)
	

	# col names for average values
	header = [ 'gain', 'power', 'kx_avg', 'ky_avg', 'kxLRO_avg', 'kyLRO_avg', 'kxvar_avg', 'kyvar_avg', 'ky_stokes' ]
	if raw:
	    header += lro_head + var_head

	# Tuple of data to be saved in the same order as header listing
	data = (self.gain, self.power, self.k_avg, self.LRO_avg, self.variance_k_avg, self.stokes_avg)

	if raw:
	    # Except first column (gain), append raw data loaded from each file
	    data += (self.LRO_raw[:,1:], self.var_raw[:,1:])

	# data variable must be reshaped so it can be stacked column ways when saved
	def my_shape(x):
	    if x.ndim > 1:
		return x
	    else:
		return reshape(x, (-1, 1))

	save_data = hstack(map(my_shape, data))

	with open(filename, 'w') as fh:
	    fh.write(comments)
	    fh.write( delimiter.join(header) + '\n')
	    savetxt( fh, save_data, fmt='%.6e', delimiter=delimiter)

	print "Data saved to file %s" % filename

def lagrange(points):
   """ Lagrange Interpolation through a set of points """
   tmp = scipy.poly1d([0])
   result=scipy.poly1d([0])
   for i in points.keys():
      numerator=scipy.poly1d([1])
      denom = 1.0
      for j in points.keys():
         if (i != j):
            tmp = scipy.poly1d([1,-j])
            numerator = numerator * tmp
            denom = denom * (i - j)
      tmp = (numerator/denom) * points.get(i)
      result = result + tmp
   return result

def main():
    data = sp.genfromtxt("web_traffic.tsv", delimiter="\t")
    plt.xkcd()

    x = data[:, 0]
    y = data[:, 1]

    x = x[~sp.isnan(y)]
    y = y[~sp.isnan(y)]

    fp1, _, _, _, _ = sp.polyfit(x, y, 1, full=True)

    # Here we try 3 degrees of freedom
    fp2, _, _, _, _ = sp.polyfit(x, y, 3, full=True)

    f1 = sp.poly1d(fp1)
    f2 = sp.poly1d(fp2)

    # We have an obvious inflection point between 3rd and 4th week
    inflection_in_hours = int(3.5 * 7 * 24)

    x_before_inflection = x[:inflection_in_hours]

    x_after_inflection = x[inflection_in_hours:]
    y_after_inflection = y[inflection_in_hours:]

    f_after = sp.poly1d(sp.polyfit(x_after_inflection, y_after_inflection, 1))

    fx = sp.linspace(0, x[-1], 1000)
    fx_after = sp.linspace(len(x_before_inflection)+1, x[-1], 1000)

    plt.scatter(x, y, s=5)
    plt.title("Web traffic over the last month.")
    plt.xlabel("Time")
    plt.ylabel("Hits/hour")
    plt.xticks([w * 7 * 24 for w in range(10)],
               ['week {}'.format(w) for w in range(10)])
    plt.autoscale(tight=True)

    plt.plot(fx, f1(fx), linewidth=2)
    plt.plot(fx, f2(fx), linewidth=2)
    plt.plot(fx_after, f_after(fx_after), linewidth=3)
    plt.legend(["d={}".format(f1.order),
                "d={}".format(f2.order),
                "d after inflection"],
               loc="upper left")
    # plt.grid(True, linestyle="-", color='0.75')
    plt.show()

def two_line(x, y):
     x1 = x[ x<3*24*7]
     x2 = x[ x>=3*24*7]
     y1 = y[ x<3*24*7]
     y2 = y[ x>=3*24*7]
     fp1 = sp.polyfit(x1, y1, 1)
     fp2 = sp.polyfit(x2, y2, 1)
     f1= sp.poly1d(fp1)
     f2 = sp.poly1d(fp2)
     plotdata(x, y)
     from matplotlib import pyplot as pl
     x_1 = sp.linspace(0, x1[-1], 1000)
     x_2 = sp.linspace(x1[-1], x2[-1], 1000)
     pl.plot(x_1, f1(x_1), linewidth=3, color="blue")
     pl.plot(x_2, f2(x_2), linewidth=3, color="red")
     pl.show()

 def fit_indices(self, indices):
     order = self.values['Order']
     if order:
         if order == -1:
             import fitfun
             def fitf(arg, y0,y1,y2,x1):
                 n=len(indices)
                 x=arg
                 ya = (y1-y0)/x1*x + y0
                 yb = (y2-y1)/(n-x1)*(x-x1)+y1
                 return ya*(x<x1)+yb*(x>=x1)
             xx = n.arange(len(indices))
             oo=fitfun.Optimizer(xx, indices)
             oo.set_function(fitf)
             oo.set_parameter_range('y0', min(indices),max(indices),0)
             oo.set_parameter_range('y1', min(indices),max(indices),0)
             oo.set_parameter_range('y2', min(indices),max(indices),0)
             oo.set_parameter_range('x1', 2.0, len(indices)-2.,len(indices)/2.)
             oo.optimize()
             #print oo.solutions
             #print 'old',indices.tolist()
             indices=fitf(arg=xx, **oo.solutions).astype('int')
             #print 'new',indices.tolist()
         else:
             #print 'old i', indices.tolist()
             x = range(len(indices))
             fit_f= poly1d( polyfit( x,indices, self.values['Order']) )
             indices = fit_f(x).round().astype('int')
             #print 'new i', indices.tolist()
     else:
         pass
     return indices

	def averaging(self, data, Start = None, End = None, N = 1, m = 3):
		'''	Усереднення між заданими вузлами.
		m	-	порядок полінома
		N	-	кількість вузлів
		'''
		x, y = data[:,0], data[:,1]
		#Обрізка в заданих межах
		if not Start is None:
			if 0 < Start < x.max():
				x, y = x[x >= Start], y[x >= Start]
		if not End is None:
			if 0 < End <= x.max():
				x, y = x[x <= End], y[x <= End]
		
		n = int(N)
		EQ = sp.poly1d( sp.polyfit(x, y, m) )
		poly_Y = EQ( x )
		Range = range(0,len(x),n)
		i = 0
		xnew, ynew = [], []
		for j in list(Range[1:])+[len(x)-1]:
			if j-i <=1:	break
			x_temp = x[i:j].mean()
			xnew.append(x_temp)
			ynew.append( (y[i:j] - poly_Y[i:j]).mean() + EQ(x_temp))
			i = j
		if self.showTmp:
			return Array(sp.array([xnew, ynew]).T, Type = data.Type, scale = data.scale),  x, y, poly_Y
		else:
			return Array(sp.array([xnew, ynew]).T, Type = data.Type, scale = data.scale)

	def AutoB_splineS(self, state, isSignal = True, senderType = 0, param = 0.95):
		'''Штучний підбір коефіцієнтів для b-сплайн інтерполяції'''
		Dict = {
			'cAutoB_splineS' : (0, 'cB_splineS',  'cB_splineStep', 'cB_splineK'),
			'sAutoB_splineS' : (1, 'sB_splineS',  'sB_splineStep', 'sB_splineK'),
			'rAutoB_splineS' : (2, 'rB_splineS',  'rB_splineStep', 'rB_splineK')}
		senderName = ''
		if isSignal:
			senderName = self.sender().objectName()
		else:
			Names = ['c', 's', 'r']
			senderName = Names[senderType]+'AutoB_splineS'
		active = (Dict[senderName][0],) + self.findChilds(QtGui.QDoubleSpinBox,Dict[senderName][1:])
		data = self.getData(active[0])
		if state:
			active[1].setEnabled(False)
			y = data[:,1]
			x = data[:,0]
			EQ = sp.poly1d( sp.polyfit(x, y, 3) )
			poly_Y = EQ( x )
			Y = y - poly_Y
			Step = float(active[2].value())
			K = float(active[3].value())

			try:
				print(str((1+Step/K**3)*param))
				active[1].setValue(sp.std(Y)**2*len(y)*(1+Step/K**2)*param)
			except:
				print("AutoB_splineS: SmoothParamError")
		else:
			active[1].setEnabled(True)

    def getPoly(self, x, y, deg):

        f = self.polyfy(x, y, deg)
        fp = sp.poly1d(f)
        error = self.error(fp, x, y)
        # logging.info(" error : %i --- function deg : %i -- function : %s ", error, deg, fp)
        return fp

def fit_linear(self, head_frac=0.15, tail_frac=0.01,
                   improvement_threshold=10, order=1, show_fit=True):
        """
        Find minimum number of reads in background set.
        Here sorted array is iteratively sliced from the end with bins having high reads values till
        sliced array shows 'good' linearity
        """
        prev_residual = None
        error_improvement = improvement_threshold+1
        sliced_sorted_vals = self.sortedArray[int(head_frac*len(self.sortedArray)):]
        rem_bins = int(0.01*len(self.sortedArray))
        while error_improvement > improvement_threshold:
            sliced_sorted_vals = sliced_sorted_vals[:-rem_bins]
            x_vals = [i for i in range(len(sliced_sorted_vals))]
            pf, residuals, rank, sv, rcond = polyfit(x_vals, sliced_sorted_vals, order, full=True)
            if prev_residual is not None:
                error_improvement = prev_residual - residuals[0]
            prev_residual = residuals[0]
        linearEq = poly1d(pf)
        x_vals = [i for i in range(len(sliced_sorted_vals))]
        pred_y_vals = linearEq(x_vals)
        if show_fit is True:
            plt.scatter(x_vals, sliced_sorted_vals, s=2, alpha=0.1, c='r', lw=0)
            plt.plot(x_vals, pred_y_vals)
            plt.show()
        self.linearMaxVal = pred_y_vals[-1]
        return True      

def norm_group(pos,trail,**kargs):
    """Takes the drain current for each file in group and builds an analysis file and works out the mean drain"""
    if "signal" in kargs:
        signal=kargs["signal"]
    else:
            signal="fluo"
    lfit=kargs["lfit"]
    rfit=kargs["rfit"]

    posfile=SA.AnalyseFile()
    posfile.metadata=pos[0].metadata
    posfile=posfile&pos[0].column(0)
    posfile.column_headers=['Energy']
    for f in pos:
        print str(f["run"])+str(f.find_col(signal))
        posfile=posfile&f.column(signal)
    posfile.add_column(lambda r:np.mean(r[1:]),"mean drain")
    ec=posfile.find_col('Energy')
    md=posfile.find_col('mean drain')
    posfile=SA.AnalyseFile(posfile)
    linearfit=scipy.poly1d(posfile.polyfit(ec,md,1,lambda x,y:lfit[0]<=x<=lfit[1]))
    posfile.add_column(lambda r:r[md]-linearfit(r[ec]),'minus linear')
    highend=posfile.mean('minus',lambda r:rfit[0]<=r[ec]<=rfit[1])
    ml=posfile.find_col('minus linear')
    posfile.add_column(lambda r:r[ml]/highend,"normalised")
    if "group_key" in kargs:
        posfile[kargs["group_key"]]=pos.key
    return posfile

    def correct_slope(self,rank=4):
        
        index = np.linspace(0,self.nsamp-1,self.nsamp)
        trends = np.zeros((self.nsamp,self.norders))
        i = 0
        '''
        for order in self.orders:
            fsubs = np.isfinite(order['flux'])
            pars = polyfit(index[fsubs][20:-20],order['flux'][fsubs][20:-20],deg=rank,w=order['uflux'][fsubs][20:-20])
            trend = poly1d(pars)
            trend_samp = trend(index)
            
            trends[:,i] = trend_samp/np.median(trend_samp)    
            i += 1

        trend_mean = np.median(trends,axis=1)

        '''
        
        fluxes = np.zeros((self.nsamp,self.norders))
        for i in np.arange(self.norders):
            flux = self.orders[i]['flux']
            flux = flux/np.median(flux)    
            fluxes[:,i] = flux
            i += 1
            
        trend_mean = np.nanmedian(fluxes,axis=1)
        fsubs = np.isfinite(trend_mean)
        pars = polyfit(index[fsubs],trend_mean[fsubs],deg=rank)
        trend_smooth = poly1d(pars)(index)
        
        for order in self.orders:
            order['flux'] /= trend_smooth