def main():
    SAMPLE_NUM = 10
    degree = 9
    x, y = sin_wgn_sample(SAMPLE_NUM)
    fig = pylab.figure(1)
    pylab.grid(True)
    pylab.xlabel('x')
    pylab.ylabel('y')
    pylab.axis([-0.1,1.1,-1.5,1.5])

    # sin(x) + noise
    # markeredgewidth mew
    # markeredgecolor mec
    # markerfacecolor mfc

    # markersize      ms
    # linewidth       lw
    # linestyle       ls
    pylab.plot(x, y,'bo',mew=2,mec='b',mfc='none',ms=8)

    # sin(x)
    x2 = linspace(0, 1, 1000)
    pylab.plot(x2,sin(2*x2*pi),'#00FF00',lw=2,label='$y = \sin(x)$')

    # polynomial fit
    reg = exp(-18)
    w = curve_poly_fit(x, y, degree,reg) #w = polyfit(x, y, 3)
    po = poly1d(w)      
    xx = linspace(0, 1, 1000)
    pylab.plot(xx, po(xx),'-r',label='$M = 9, \ln\lambda = -18$',lw=2)
    
    pylab.legend()
    pylab.show()
    fig.savefig("poly_fit9_10_reg.pdf")

def display_coeff(data=None):
    betaAll,betaErrAll, R2adjAll = measure_stamp_coeff(data = data, zernike_max_order=20)
    ind = np.arange(len(betaAll[0]))
    momname = ('M20','M22.Real','M22.imag','M31.real','M31.imag','M33.real','M33.imag')
    fmtarr = ['bo-','ro-','go-','co-','mo-','yo-','ko-']
    pl.figure(figsize=(17,13))
    for i in range(7):
        pl.subplot(7,1,i+1)
        pl.errorbar(ind,betaAll[i],yerr = betaErrAll[i],fmt=fmtarr[i])
        pl.grid()
        pl.xlim(-1,21)
        if i ==0:
            pl.ylim(-10,65)
        elif i ==1:
            pl.ylim(-5,6)
        elif i ==2:
            pl.ylim(-5,6)
        elif i == 3:
            pl.ylim(-0.1,0.1)
        elif i == 4:
            pl.ylim(-0.1,0.1)
        elif i ==5:
            pl.ylim(-100,100)
        elif i == 6:
            pl.ylim(-100,100)
        pl.xticks(ind,('','','','','','','','','','','','','','','','','','','',''))
        pl.ylabel(momname[i])
    pl.xticks(ind,('0','1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19'))
    pl.xlabel('Zernike Coefficients')
    return '--- done ! ----'

def drawPr(tp,fp,tot,show=True):
    """
        draw the precision recall curve
    """
    det=numpy.array(sorted(tp+fp))
    atp=numpy.array(tp)
    afp=numpy.array(fp)
    #pylab.figure()
    #pylab.clf()
    rc=numpy.zeros(len(det))
    pr=numpy.zeros(len(det))
    #prc=0
    #ppr=1
    for i,p in enumerate(det):
        pr[i]=float(numpy.sum(atp>=p))/numpy.sum(det>=p)
        rc[i]=float(numpy.sum(atp>=p))/tot
        #print pr,rc,p
    ap=0
    for c in numpy.linspace(0,1,num=11):
        if len(pr[rc>=c])>0:
            p=numpy.max(pr[rc>=c])
        else:
            p=0
        ap=ap+p/11
    if show:
        pylab.plot(rc,pr,'-g')
        pylab.title("AP=%.3f"%(ap))
        pylab.xlabel("Recall")
        pylab.ylabel("Precision")
        pylab.grid()
        pylab.show()
        pylab.draw()
    return rc,pr,ap

def plot_sphere_x( s, fname ):
  """ put plot of ionization fractions from sphere `s` into fname """

  plt.figure()
  s.Edges.units = 'kpc'
  s.r_c.units = 'kpc'
  xx = s.r_c
  L = s.Edges[-1]

  plt.plot( xx, np.log10( s.xHe1 ),
            color='green', ls='-', label = r'$x_{\rm HeI}$' )
  plt.plot( xx, np.log10( s.xHe2 ),
            color='green', ls='--', label = r'$x_{\rm HeII}$' )
  plt.plot( xx, np.log10( s.xHe3 ),
            color='green', ls=':', label = r'$x_{\rm HeIII}$' )

  plt.plot( xx, np.log10( s.xH1 ),
            color='red', ls='-', label = r'$x_{\rm HI}$' )
  plt.plot( xx, np.log10( s.xH2 ),
            color='red', ls='--', label = r'$x_{\rm HII}$' )

  plt.xlim( -L/20, L+L/20 )
  plt.xlabel( 'r_c [kpc]' )

  plt.ylim( -4.5, 0.2 )
  plt.ylabel( 'log 10 ( x )' )

  plt.grid()
  plt.legend(loc='best', ncol=2)
  plt.tight_layout()
  plt.savefig( 'doc/img/x_' + fname )

    def plot_number_alteration_by_tissue(self, fontsize=10, width=0.9):
        """Plot number of alterations

        .. plot::
            :width: 100%
            :include-source:

            from gdsctools import *
            data = gdsctools_data("test_omnibem_genomic_alterations.csv.gz")
            bem = OmniBEMBuilder(data)
            bem.filter_by_gene_list(gdsctools_data("test_omnibem_genes.txt"))
            bem.plot_number_alteration_by_tissue()

        """
        count = self.unified.groupby(['TISSUE_TYPE'])['GENE'].count()
        try:
            count.sort_values(inplace=True, ascending=False)
        except:
            count.sort(inplace=True, ascending=False)
        count.plot(kind="bar", width=width)
        pylab.grid()
        pylab.xlabel("Tissue Type", fontsize=fontsize)
        pylab.ylabel("Total number of alterations in cell lines",
                     fontsize=fontsize)
        try:pylab.tight_layout()
        except:pass

    def plot_roc(self, roc=None):
        """Plot ROC curves

        .. plot::
            :include-source:
            :width: 80%

            from dreamtools import rocs
            r = rocs.ROC()
            r.scores = [.9,.5,.6,.7,.1,.2,.6,.4,.7,.9, .2]
            r.classes = [1,0,1,0,0,1,1,0,0,1,1]
            r.plot_roc()

        """
        if roc == None:
            roc = self.get_roc()
        from pylab import plot, xlim, ylim ,grid, title, xlabel, ylabel
        x = roc['fpr']
        plot(x, roc['tpr'], '-o')
        plot([0,1], [0,1],'r')
        ylim([0, 1])
        xlim([0, 1])
        grid(True)
        title("ROC curve (AUC=%s)" % self.compute_auc(roc))
        xlabel("FPR")
        ylabel("TPR")

def plot_signal(x,y,title,labelx,labley,position):
    pylab.subplot (9, 1, position)
    pylab.plot(x,y)
    pylab.title(title)
    pylab.xlabel(labelx)
    pylab.ylabel(labley)
    pylab.grid(True)

	def update(self):
		if self.pose != []:
			plt.figure(1)
			clf()
			self.fig1 = plt.figure(num=1, figsize=(self.window_size, \
				self.window_size), dpi=80, facecolor='w', edgecolor='w')
			title (self.title)			
			xlabel('Easting [m]')
			ylabel('Northing [m]')
			axis('equal')
			grid (True)
			poseT = zip(*self.pose)	
			pose_plt = plot(poseT[1],poseT[2],'#ff0000')

			if self.wptnav != []:
				mode = self.wptnav[-1][MODE]

				if not (self.wptnav[-1][B_E] == 0 and self.wptnav[-1][B_N] == 0 and self.wptnav[-1][A_E] == 0 and self.wptnav[-1][A_N] == 0):
					b_dot = plot(self.wptnav[-1][B_E],self.wptnav[-1][B_N],'ro',markersize=8)
					a_dot = plot(self.wptnav[-1][A_E],self.wptnav[-1][A_N],'go',markersize=8)
					ab_line = plot([self.wptnav[-1][B_E],self.wptnav[-1][A_E]],[self.wptnav[-1][B_N],self.wptnav[-1][A_N]],'g')
					target_dot = plot(self.wptnav[-1][TARGET_E],self.wptnav[-1][TARGET_N],'ro',markersize=5)

				if mode == -1:
					pose_dot = plot(self.wptnav[-1][POSE_E],self.wptnav[-1][POSE_N],'b^',markersize=8)
				elif mode == 1:
					pose_dot = plot(self.wptnav[-1][POSE_E],self.wptnav[-1][POSE_N],'bs',markersize=8)
				elif mode == 2:
					pose_dot = plot(self.wptnav[-1][POSE_E],self.wptnav[-1][POSE_N],'bo',markersize=8)

		if self.save_images:
			self.fig1.savefig ('plot_map%05d.jpg' % self.image_count)
			self.image_count += 1
		draw()

def plot_spectrum():
    #get the data...    
    a_0=struct.unpack('>1024l',fpga.read('even',1024*4,0))
    a_1=struct.unpack('>1024l',fpga.read('odd',1024*4,0))

    interleave_a=[]

    for i in range(1024):
        interleave_a.append(a_0[i])
        interleave_a.append(a_1[i])

    pylab.figure(num=1,figsize=(10,10))
    pylab.ioff()
    pylab.plot(interleave_a)
    #pylab.semilogy(interleave_a)
    pylab.title('Integration number %i.'%prev_integration)
    pylab.ylabel('Power (arbitrary units)')
    pylab.grid()
    pylab.xlabel('Channel')
    pylab.xlim(0,2048)
    pylab.ioff()

    pylab.hold(False)
    pylab.show()
    pylab.draw()

def pr_curve(mean_precision, mean_recall, title='20-newsgroups', prefix='', suffix='', save_path=''):
    """
    Parameters
    ----------

    mean_precision: numpy.ndarry，一维
        记录平均精度的向量.
    mean_recall: numpy.ndarray，一维
        记录平均召回率向量.
    prefix, duffix: 同rsm_query()
    save_path: str, 保存图片的目录，注意一定要使用绝对路径，支持以'/'结尾或没有
    """  
    #color_list = ['blue', 'red', 'green', 'cyan', 'yellow', 'black', 'magenta', (0.5,0.5,0.5)]
    y_vector = mean_precision * 100.
    x_vector = mean_recall * 100 #** (1./6.)

    pylab.figure(figsize=(8, 8))
    pylab.grid() #在做标系中显示网格
    pylab.plot(x_vector, y_vector, label='$r-p curve$', color='blue', linewidth=1)

    pylab.xlabel('recall(%)')
    pylab.ylabel('precision(%)')
    pylab.title(title)
    #pylab.xlim(0., 60) #x轴长度限制
    #pylab.ylim(0., 30) #y轴长度限制
    pylab.legend() #在图像中显示标记说明
    #pylab.show() # 显示图像
    if save_path != '':
        if platform.system() == 'Linux' and save_path[-1] != '/':
            save_path += '/'
        if platform.system() == 'Windows' and save_path[-1] != '\\':
            save_path += '\\'
    save_path += 'plots/'
    check_path(save_path) #检测路径，如果目录不存在，支持递归建立
    pylab.savefig(save_path + prefix + 'r-p' + suffix + '.png', dpi=240) #保存图像，可以人为指定所保存的图像的分辨率

	def analyze(self,igroup):
		# to get the sourcedata
		database=self.gettmpdata('database');
		spectra=database[0]['resultdatatablegroups'][igroup];
		spectranew=spectra.getemptyinstance();
		
		#the real calculation
		spectraall=database[0]['datalist'][0];
		spectraall.plot();
		#spectraall.show();
		pylab.grid(True);

		formulastr=eval(self['formulastr'].get());
		
		for k in spectra.keys():
			spect=spectra[k];
			spectnew=spect.copyxy();
			spectnew.update(formulastr);
			spectnew.log({"Operation":"set column value","formulastr":formulastr});
			spectranew.insert(spectnew,k);
					
		spectra=database[0]['resultdatatablegroups'][igroup]=spectranew;
		#pylab.figure();
		XpyFigure();
		spectraall.plot();
		spectranew.plot('o');

	def test_seq(self):
		sleep(1)
		
		self.cmd_logger(0)
		self.takeoffpub.publish(Empty())	;print 'takeoff' #takeoff
		sleep(12); print '4'; sleep(1); print '3'; sleep(1); print '2'; sleep(1); print '1'; sleep(1)
		
		self.cmd_logger(0)
		self.twist.linear.z = self.vzcmd
		self.cmd_logger(self.vzcmd)
		self.cmdpub.publish(self.twist)		;print 'vzcmd' #set vzcmd
		sleep(self.vzdur)			#wait for vzdur
		
		self.cmd_logger(self.vzcmd)
		self.clear_twist()
		self.cmd_logger(0)
		self.cmdpub.publish(self.twist)		;print 'clear vz' #clear vz
		sleep(4)
		
		self.cmd_logger(0)
		self.landpub.publish(Empty())		;print 'land' #land
		sleep(1)
		
		if not raw_input('show and save?') == 'n':
			pl.xlabel('time (s)')
			pl.ylim(0,2400)
			pl.plot(self.cmd_log['tm'],self.cmd_log['cmd'], 'b-s')
			pl.plot(self.nd_log['tm'],self.nd_log['vz'], 'g-+')
			pl.plot(self.nd_log['tm'],self.nd_log['al'], 'r-+')
			pl.grid(True)
			pl.show()

def hist_Ne(sel_coinc_ids, coords, data, n):
    global histNe2
    c_index = data.root.coincidences.c_index
    observ = data.root.coincidences.observables
    core_rec = data.root.core_reconstructions.reconstructions
    
    histNe = core_rec.readCoordinates(coords, field='reconstructed_shower_size')
    #histNe = [x for x in histNe if x > 0]  # for showersize smaller than 0 
    
    d = 10**10.2
    
    histNe2 = [x*d for x in histNe]
    #histNe *= 10 ** 10.2
    
      
    pylab.hist(np.log10(histNe2), 100, log=True) # histtype="step"
       
    pylab.xlabel('Showerenergy log(eV)')
    pylab.ylabel('count')
    pylab.title('showersize bij N==%s' %n)
    pylab.ylim(ymin=1)
    pylab.grid(True)
    pylab.show()
    
    return histNe2

def plot_1d(profile, title_fig, title_x, title_y):
    # pylab.figure()
    pylab.plot(profile[0], profile[1])
    pylab.xlabel(title_x)
    pylab.ylabel(title_y)
    pylab.title(title_fig)
    pylab.grid(True)

    def __call__(self, n):
        if len(self.f.shape) == 3:
            # f = f[x,v,t], 2 dim in phase space
            ft = self.f[n,:,:]
            pylab.pcolormesh(self.X, self.V, ft.T, cmap = 'jet')
            pylab.colorbar()
            pylab.clim(0,0.38) # for Landau test case
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$v$', fontsize = 18)
            pylab.title('$N_x$ = %d, $N_v$ = %d, $t$ = %2.1f' % (self.x.N, self.v.N, self.it*self.t.width))
            pylab.savefig(self.path + self.filename)
            pylab.clf()
            return None

        if len(self.f.shape) == 2:
            # f = f[x], 1 dim in phase space
            ft = self.f[n,:]
            pylab.plot(self.x.gridvalues,ft,'ob')
            pylab.grid()
            pylab.axis([self.xmin, self.xmax, self.ymin, self.ymax])
            pylab.xlabel('$x$', fontsize = 18)
            pylab.ylabel('$f(x)$', fontsize = 18)
            pylab.savefig(self.path + self.filename)
            return None

def drunkTest(numTrials = 1000):
    #stepsTaken = [10, 100, 1000, 10000]
    stepsTaken = 1000
    
    for dClass in (UsualDrunk, ColdDrunk, EDrunk, PhotoDrunk, DDrunk):
        #initialize field
        field = Field()
        origin = Location(0, 0)
        
        # initialize drunk 
        drunk = dClass('Drunk')
        field.addDrunk(drunk, origin)

        x_pos, y_pos = [], [] # initialize to empty
        x, y = 0.0, 0.0
        
        for trial in range(numTrials): # trials 
            x, y = walkVector(field, drunk, stepsTaken)
            x_pos.append(x)
            y_pos.append(y)
            
        #pylab.plot(x_pos, y_pos, 'ro', s=5,
        #           label = dClass.__name__)
        pylab.scatter(x_pos, y_pos,s=5, color='red')
        pylab.title(str(dClass))
        pylab.xlabel('x')
        pylab.grid()
        pylab.xlim(-100, 100)
        pylab.ylim(-100,100)
        pylab.ylabel('y')
        pylab.show()

def plot_ch():
    for job in jobs_orig:
        print "plane of", job.path
        pylab.clf()
        x_center = int((job.position(0)[0] + job.position(1)[0])/2)
        x_final = 50 + x_center
        #plane = np.concatenate((job.plane(y=50)[:, x_final:], 
        #                        job.plane(y=50)[:, :x_final]), axis=1)
        plane = job.plane(y=50)
        myplane = plane[plane < 0.0]
        p0 = myplane.min()
        p12 = np.median(myplane)
        p14 = np.median(myplane[myplane<p12])
        p34 = np.median(myplane[myplane>p12])
        p1 = myplane.max()
        contour_values = (p0, p14, p12, p34, p1)
        pylab.title(r'$u_x=%.4f,\  D_{-}=%.4f,\  D_{+}=%.4f,\ ch=%i$ ' %
                    (job.u_x, job.D_minus, job.D_plus, job.ch_objects))
        car = pylab.imshow(plane, vmin=-0.001, vmax=0.0, 
                           interpolation='nearest')
        pylab.contour(plane, contour_values, linestyles='dashed', 
                                             colors='white')
        pylab.grid(True)
        pylab.colorbar(car)
        #imgfilename = 'plane_r20-y50-u_x%.4fD%.4fch%03i.png' % \
        #              (job.u_x, job.D_minus, job.ch_objects)
        imgfilename = 'plane_%s.png' % job.job_id
        pylab.savefig(imgfilename)

def main():
    # The Wiener process parameter.
    delta = 2
    # Total time.
    T = 10.0
    # Number of steps.
    N = 1000
    # Time step size
    dt = T/N
    # Initial value of x.
    x0 = 100.0

    time0 = time.time()
    xa = random_walk(x0, N, dt, delta)
    time1 = time.time()
    xb = random_walk_loop(x0, N, dt, delta)
    time2 = time.time()

    print "numpy cumsum: %.3g seconds" % (time1-time0)
    print "python loop:  %.3g seconds" % (time2-time1)
    print "ratio: %.2f" % ((time2-time1)/(time1-time0))

    t = numpy.linspace(0.0, N*dt, N+1)
    plot(t, xa, 'g')
    plot(t, xb, 'b')
    xlabel('t', fontsize=16)
    ylabel('x',fontsize=16)
    grid(True)
    show()

def TestOverAlpha():
    nDim = 5
    numOfParticles = 10
    maxIteration = 2000
    minX = array([-100.0]*nDim)
    maxX = array([100.0]*nDim)
    maxV = 1.0*(maxX - minX)
    minV = -1.0*maxV
    numOfTrial = 10
    intDim = 4
    alpha = 0.3
    while alpha<1.0:
        gBest = array([0.0]*maxIteration)
        for i in xrange(numOfTrial):
            p1 = AUPSO.PSOProblem(nDim, numOfParticles, maxIteration, minX, maxX, minV, maxV, AUPSO.Sphere,intDim,alpha)
            p1.run()
            gBest = gBest + p1.gBestArray[:maxIteration]
        gBest = gBest / numOfTrial
        pylab.plot(range(maxIteration), gBest,label='alpha='+str(alpha))
        print 'alpha = ', alpha
        alpha += 0.3
    print 'now drawing'
    pylab.title('$G_{best}$ over 20 trials'+' intDim='+str(intDim))
    pylab.xlabel('The $N^{th}$ Iteratioin')
    pylab.ylabel('Average gBest over '+str(numOfTrial)+' runs')
    pylab.grid(True)
    pylab.yscale('log')
    ylim = [-6, 1]
    ystep = 1.0
#    pylab.ylim(ylim[0], ylim[1])
#    yticks = linspace(ylim[0], ylim[1], int((ylim[1]-ylim[0])/ystep+1))
#    pylab.yticks(tuple(yticks), tuple(map(str,yticks)))
    pylab.legend(loc='lower left')
    pylab.show()

 def testTelescope(self):
     import matplotlib
     matplotlib.use('AGG')
     import matplotlib.mlab as ml
     import pylab as pl
     import time        
     w0 = 8.0
     k = 2*np.pi/3.0
     gb = GaussianBeam(w0, k)
     lens = ThinLens(150, 150)
     gb2 = lens*gb
     self.assertAlmostEqual(gb2._z0, gb._z0 + 2*150.0)
     lens2 = ThinLens(300, 600)
     gb3 = lens2*gb2
     self.assertAlmostEqual(gb3._z0, gb2._z0 + 2*300.0)
     self.assertAlmostEqual(gb._w0, gb3._w0/2.0)
     z = np.arange(0, 150)
     z2 = np.arange(150, 600)
     z3 = np.arange(600, 900)
     pl.plot(z, gb.w(z, k), z2, gb2.w(z2, k), z3, gb3.w(z3, k))
     pl.grid()
     pl.xlabel('z')
     pl.ylabel('w')
     pl.savefig('testTelescope1.png')
     time.sleep(0.1)
     pl.close('all')