Python scipy.poly1d函数代码示例

    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)