def PSTH(self):
	
			
		TimeRes = np.array([0.1,0.25,0.5,1,2.5,5.0,10.0,25.0,50.0,100.0])

		Projection_PSTH = np.zeros((2,len(TimeRes)))
		for i in range(0,len(TimeRes)):
			Data_Hist,STA_Hist,Model_Hist,B = Hist(TimeRes[i])
			data = Data_Hist/np.linalg.norm(Data_Hist)
			sta = STA_Hist/np.linalg.norm(STA_Hist)
			model = Model_Hist/np.linalg.norm(Model_Hist)
			Projection_PSTH[0,i] = np.dot(data,sta)
			Projection_PSTH[1,i] = np.dot(data,model)
			
		import matplotlib.font_manager as fm
		
		plt.figure()
		plt.semilogx(TimeRes,Projection_PSTH[0,:],'gray',TimeRes,Projection_PSTH[1,:],'k',
			     linewidth=3, marker='o', markersize = 12)
		plt.xlabel('Time Resolution, ms',fontsize=25)
		plt.xticks(fontsize=25)
		#plt.axis["right"].set_visible(False)
		plt.ylabel('Projection onto PSTH',fontsize=25)
		plt.yticks(fontsize=25)
		prop = fm.FontProperties(size=20)
		plt.legend(('1D model','2D model'),loc='upper left',prop=prop)
		plt.tight_layout()
		plt.show()

def make_corr1d_fig(dosave=False):
    corr = make_corr_both_hemi()
    lw=2; fs=16
    pl.figure(1)#, figsize=(8, 7))
    pl.clf()
    pl.xlim(4,300)
    pl.ylim(-400,+500)    
    lambda_titles = [r'$20 < \lambda < 30$',
                     r'$30 < \lambda < 40$',
                     r'$\lambda > 40$']
    colors = ['blue','green','red']
    for i in range(3):
        corr1d, rcen = corr_1d_from_2d(corr[i])
        ipdb.set_trace()
        pl.semilogx(rcen, corr1d*rcen**2, lw=lw, color=colors[i])
        #pl.semilogx(rcen, corr1d*rcen**2, 'o', lw=lw, color=colors[i])
    pl.xlabel(r'$s (Mpc)$',fontsize=fs)
    pl.ylabel(r'$s^2 \xi_0(s)$', fontsize=fs)    
    pl.legend(lambda_titles, 'lower left', fontsize=fs+3)
    pl.plot([.1,10000],[0,0],'k--')
    s_bao = 149.28
    pl.plot([s_bao, s_bao],[-9e9,+9e9],'k--')
    pl.text(s_bao*1.03, 420, 'BAO scale')
    pl.text(s_bao*1.03, 370, '%0.1f Mpc'%s_bao)
    if dosave: pl.savefig('xi1d_3bin.pdf')

def analyze(title, x, y, func, func_title):
    print('-' * 80)
    print(title)
    print('x: %s:%s %s' % (list(x.shape), x.dtype, [x.min(), x.max()]))
    print('y: %s:%s %s' % (list(y.shape), y.dtype, [y.min(), y.max()]))

    popt, pcov = curve_fit(func, x, y)
    print('popt=%s' % popt)
    print('pcov=\n%s' % pcov)

    a, b = popt
    print('a=%e' % a)
    print('b=%e' % b)
    print(func_title(a, b))
    xf = np.linspace(x.min(), x.max(), 100)
    yf = func(xf, a, b)
    print('xf: %s:%s %s' % (list(xf.shape), xf.dtype, [xf.min(), xf.max()]))
    print('yf: %s:%s %s' % (list(yf.shape), yf.dtype, [yf.min(), yf.max()]))

    plt.title(func_title(a, b))
    # plt.xlim(0, x.max())
    # plt.ylim(0, y.max())
    plt.semilogx(x, y, label='data')
    plt.semilogx(xf, yf, label='fit')
    plt.legend(loc='best')
    plt.savefig('%s.png' % title)
    plt.close()

def pore_size_distribution(network, fig=None):
    r"""
    Plot the pore and throat size distribution which is the accumulated
    volume vs. the diameter in a semilog plot

    Parameters
    ----------
    network : OpenPNM Network object

    """
    if fig is None:
        fig = _plt.figure()
    dp = network['pore.diameter']
    Vp = network['pore.volume']
    dt = network['throat.diameter']
    Vt = network['throat.volume']
    dmax = max(max(dp), max(dt))
    steps = _sp.linspace(0, dmax, 100, endpoint=True)
    vals = _sp.zeros_like(steps)
    for i in range(0, len(steps)-1):
        temp1 = dp > steps[i]
        temp2 = dt > steps[i]
        vals[i] = sum(Vp[temp1]) + sum(Vt[temp2])
    yaxis = vals
    xaxis = steps
    _plt.semilogx(xaxis, yaxis, 'b.-')
    _plt.xlabel('Pore & Throat Diameter (m)')
    _plt.ylabel('Cumulative Volume (m^3)')
    return fig

    def check_models(self):
        plt.figure('Bandgap narrowing')
        Na = np.logspace(12, 20)
        Nd = 0.
        dn = 1e14
        temp = 300.

        for author in self.available_models():
            BGN = self.update(Na=Na, Nd=Nd, nxc=dn,
                              author=author,
                              temp=temp)

            if not np.all(BGN == 0):
                plt.plot(Na, BGN, label=author)

        test_file = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            'Si', 'check data', 'Bgn.csv')

        data = np.genfromtxt(test_file, delimiter=',', names=True)

        for name in data.dtype.names[1:]:
            plt.plot(
                data['N'], data[name], 'r--',
                label='PV-lighthouse\'s: ' + name)

        plt.semilogx()
        plt.xlabel('Doping (cm$^{-3}$)')
        plt.ylabel('Bandgap narrowing (K)')

        plt.legend(loc=0)

def experiment_plot( ctr, trials, success ):
	"""
	Pass in the ctr, trials and success returned
	by the `experiment` function and plot
	the Cumulative Number of Turns For Each Arm and
	the CTR's Convergence Plot side by side
	"""
	T, K = trials.shape
	n = np.arange(T) + 1
	fig = plt.figure( figsize = ( 14, 7 ) )

	plt.subplot(121)	
	for i in range(K):
		plt.loglog( n, trials[ :, i ], label = "arm {}".format(i + 1) )

	plt.legend( loc = "upper left" )
	plt.xlabel("Number of turns")
	plt.ylabel("Number of turns/arm")
	plt.title("Cumulative Number of Turns For Each Arm")

	plt.subplot(122)
	for i in range(K):
		plt.semilogx( n, np.zeros(T) + ctr[i], label = "arm {}'s CTR".format( i + 1 ) )

	plt.semilogx( n, ( success[ :, 0 ] + success[ :, 1 ] ) / n, label = "CTR at turn t" )

	plt.axis([ 0, T, 0, 1 ] )
	plt.legend( loc = "upper left" )
	plt.xlabel("Number of turns")
	plt.ylabel("CTR")
	plt.title("CTR's Convergence Plot")

	return fig

def plot(x, y, semilogx = False):
	if semilogx:
		pylab.semilogx(x, y)
	else:
		pylab.plot(x, y)
		
	pylab.grid(True)
	pylab.show()

def deg_coll_plots(dirlist=None, **kwargs):
    if kwargs.get("perzeta",False)==True:
        zetalist = []
        for dirnames in dirlist:
          zetalist.append(zetaper(w_dir=dirnames,nrad=kwargs.get('nrad',638)))
          #f1 = plt.figure()
        plt.xlabel(r'Distance Along Field Line : s [AU]',labelpad=6)
        plt.ylabel(r'$\Delta \zeta [\%]$',labelpad=10)
        plt.plot(zetalist[0][0],zetalist[0][1],'k-.')
        plt.plot(zetalist[1][0],zetalist[1][1],'k--')
        plt.plot(zetalist[2][0],zetalist[2][1],'k')
        plt.axis([0.0,150.0,0.0,40.0])
        plt.minorticks_on()
    else:
        magalf = np.linspace(19.67,19.67,100)
        magfast = np.linspace(12.53,12.53,100)
        
        Msarr = [20.0,25.0,30.0,45.0,50.0,60.0]
        alfMs = [20.44,23.27,26.05, 26.84, 28.86,31.45]
        fastMs = [15.57,19.94,23.40,24.35,26.03,29.25]
        asymMs=[]
        Alparr = [0.55,0.60,0.65]
        alfAlp =[31.35,23.73,20.93]
        fastAlp=[29.25,21.48,17.28]
        asymAlp =[]
        Denarr =[3.0e-14,5.0e-14,1.0e-13,5.0e-13]
        alfDen =[30.30,26.05,22.91,20.42]
        fastDen =[27.79,23.40,19.47,15.03]
        asymDen=[]
        Betarr=[1.0,3.0,5.0]
        alfBet =[21.98,25.63,26.05]
        fastBet=[16.45,20.68,23.40]
        asymBet=[]

        #nice_plots()
        #f1 = plt.figure()
        #ax1 = f1.add_subplot(111)
        
        plt.ylabel(r"Opening Angle : $\phi^\circ$",labelpad=10)
        plt.minorticks_on()
        if kwargs.get('cplot',False)==True:
            dum1 = np.linspace(1.0e-14,1.0e-12,100)
            ms1 = plt.semilogx(Denarr,alfDen,'r*')
            ms2 = plt.semilogx(Denarr,fastDen,'bo')
            plt.xlabel(r"Base Density : $\rho_0$ [g/cm$^{3}$]")
            plt.semilogx(dum1,magalf,'r--',dum1,magfast,'b--')
            plt.axis([1.0e-14,1.0e-12,10.0,35.0])
           # plt.legend([ms1,ms2],[r'Alfven point',r'Fast point'],loc='lower left')
        else:
            dum1 = np.linspace(10.0,70.0,100)
            ms1 = plt.plot(Msarr,alfMs,'k*')
            ms2 = plt.plot(Msarr,fastMs,'ko')
            plt.xlabel(r"Stellar Mass : $M_{*}$ [$M_{\odot}$]",labelpad=6)
            plt.plot(dum1,magalf,'k--',dum1,magfast,'k')
            plt.axis([10.0,70.0,10.0,35.0])

def makePlot(xVals, yVals, title, xLabel, yLabel, style, logX=False, logY=False):
    """用给定的标题和标签绘制xVals和yVals
    """
    plt.figure()
    plt.title(title)
    plt.xlabel(xLabel)
    plt.ylabel(yLabel)
    plt.plot(xVals, yVals, style)
    if logX:
        plt.semilogx()
    if logY:
        plt.semilogy()

def plot_ac_to_compare(filter_data_from_spice, coef_b, coef_a, FS, visualise=False):
    if visualise:
        # get and plot data from spice
        w1, h1 = parse_data(filter_data_from_spice)
        # get and plot data from coesfs
        w, h = signal.freqz(coef_b, coef_a)
        plt.figure()
        plt.plot(w1, h1)
        plt.semilogx(w / max(w) * FS / 2, 20 * np.log10(abs(h)))
        plt.xlabel("Frequency [ Hz ]")
        plt.ylabel("Amplitude [dB]")
        plt.margins(0, 0.1)
        plt.grid(which="both", axis="both")
        plt.show()

def showErrorBars(minExp, maxRxp, numTrials):
    """绘制误差图
    """
    means, sds = [], []
    xVals = []
    for exp in range(minExp, maxRxp+1):
        xVals.append(2**exp)
    for numFlips in xVals:
        fracHeads, mean, sd = flipSim(numFlips, numTrials)
        means.append(mean)
        sds.append(sd)
    plt.errorbar(xVals, means, yerr=2*np.array(sds))
    plt.semilogx()
    plt.title('Mean Fraction of Heads (' + str(numTrials) + ' Trials)')
    plt.xlabel('Number of flips per trial')
    plt.ylabel("Fraction of heads & 95% confidence")

def plot_result(fname='results/experiment_run.npy'):
    """ Plot beta/density overview
    """
    data = np.load(fname)

    plt.semilogx(*zip(*data), marker='o', ls='')

    plt.title(r'Influence of $\beta$ on spiral-tip density')
    plt.xlabel(r'$\beta$')
    plt.ylabel(r'spiral-tip density')

    img_dir = 'images'
    if not os.path.isdir(img_dir):
        os.makedirs(img_dir)
    plt.savefig(
        os.path.join(img_dir, 'beta_overview.png'),
        bbox_inches='tight', dpi=300)

def check_klaassen():
    '''compares to values taken from www.PVlighthouse.com.au'''
    a = Mobility('Si')
    a.change_model('klaassen1992')

    print('''The model disagrees at low tempeature owing to dopant\
           ionisation\
           I am unsure if mobility should take ionisated dopants or\
            non ionisaed\
           most likley it should take both, currently it only takes one''')

    dn = np.logspace(10, 20)
    # dn = np.array([1e14])
    Nd = 1e14
    Na = 0

    folder = os.path.join(
        os.path.dirname(__file__), 'Si', 'test_mobility_files')
    fnames = [r'Klassen_1e14_dopants.dat',
              r'Klassen_1e14_temp-450.dat']
    print(os.path.isdir(folder))

    for temp, f_name in zip([300, 450], fnames):

        plt.figure('Mobility - Klaassen: Deltan at ' + str(temp))

        plt.plot(dn, a.hole_mobility(dn, Na, Nd, temp=temp),
                 'r-',
                 label='hole-here')
        plt.plot(dn, a.electron_mobility(dn, Na, Nd, temp=temp),
                 'b-',
                 label='electron-here')

        print(f_name)
        data = np.genfromtxt(os.path.join(folder, f_name), names=True)

        plt.plot(data['deltan'], data['uh'], 'b--',
                 label='hole - PV-lighthouse')
        plt.plot(data['deltan'], data['ue'], 'r--',
                 label='electron - PV-lighthouse')
        plt.legend(loc=0, title='Mobility from')

        plt.semilogx()
        plt.xlabel(r'$\Delta$n (cm$^{-3}$)')
        plt.xlabel(r'Moblity  (cm$^2$V$^{-1}$s$^{-1}$)')

def cdf(x,color='b',label=" ",lw=1,xlabel="x",ylabel="CDF",logx=False):
    """ plot the cumulative density function of x

    Parameters
    ----------

    x :  np.array  (N)
    color : string
        color symbol
    label : string
        label
    lw: float
        linewidth
    xlabel : string
        xlabel
    ylabel : string
        ylabel

    Examples
    --------

    .. plot::
        :include-source:

        >>> from matplotlib.pyplot import *
        >>> import pylayers.util.pyutil as pyu
        >>> from scipy import *
        >>> import matplotlib.pylab as plt
        >>> x = randn(100)
        >>> pyu.cdf(x)

    """
    x  = np.sort(x)
    n  = len(x)
    x2 = np.repeat(x, 2)
    y2 = np.hstack([0.0, np.repeat(np.arange(1,n) / float(n), 2), 1.0])
    if logx:
        plt.semilogx(x2,y2,color=color,label=label,linewidth=lw)
    else:
        plt.plot(x2,y2,color=color,label=label,linewidth=lw)
    plt.legend()
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)

def check_dorkel():
    '''compares to values taken from www.PVlighthouse.com.au'''

    a = Mobility('Si')
    a.change_model(author='dorkel1981')

    dn = np.logspace(10, 20)
    # dn = np.array([1e14])
    Nd = 1e14
    Na = 0

    folder = os.path.join(
        os.path.dirname(__file__), 'Si', 'test_mobility_files')
    # file name and temp its at
    compare = [
        ['dorkel_1e14_carriers.dat', 300],
        ['dorkel_1e14_temp-450.dat', 450],
    ]

    for comp in compare:

        plt.figure('Mobility - Dorkel: Deltan at ' + str(comp[1]))

        plt.plot(dn, a.hole_mobility(dn, Na, Nd, temp=comp[1]),
                 'r-',
                 label='hole-here')
        plt.plot(dn, a.electron_mobility(dn, Na, Nd, temp=comp[1]),
                 'b-',
                 label='electron-here')

        data = np.genfromtxt(os.path.join(folder, comp[0]), names=True)

        plt.plot(data['deltan'], data['uh'], 'b--',
                 label='hole - PV-lighthouse')
        plt.plot(data['deltan'], data['ue'], 'r--',
                 label='electron - PV-lighthouse')
        plt.legend(loc=0, title='Mobility from')

        plt.semilogx()
        plt.xlabel(r'$\Delta$n (cm$^{-3}$)')
        plt.xlabel(r'Moblity  (cm$^2$V$^{-1}$s$^{-1}$)')

    def check_models(self):
        '''
        Plots a check of the modeled data against Digitised data from either
        papers or from other implementations of the model.
        '''
        plt.figure('Ionised impurities')

        iN_imp = N_imp = np.logspace(15, 20)

        dn = 1e10
        temp = 300

        for impurity in ['phosphorous', 'arsenic']:

            iN_imp = self.update_dopant_ionisation(N_imp,
                                                   dn,
                                                   impurity,
                                                   temp, author=None)

            if not np.all(iN_imp == 0):
                plt.plot(
                    N_imp, iN_imp / N_imp * 100, label='Altermatt: ' + impurity)

        test_file = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            'Si', 'check data', 'donors.csv')

        data = np.genfromtxt(test_file, delimiter=',', skip_header=1)

        for i in range(0, (data.shape[1] + 2) / 2, 2):
            # print i
            plt.plot(
                data[:, i], data[:, i + 1] * 100, 'r.',
                label='Digitised data')

        plt.semilogx()
        plt.xlabel('Impurity (cm$^{-3}$)')
        plt.ylabel('Fraction of Ionised impurities (%)')

        plt.legend(loc=0)

def expt1():
	"""
	Experiment 1: Chooses the result files and generates figures
	"""
	# filename = sys.argv[1]
	result_file = "./expt1.txt"
	input_threads, input_sizes, throughputs, resp_times \
		= parse_output(result_file) 

	throughputs_MiB = [tp/2**20 for tp in throughputs]

	fig1 = pl.figure()
	fig1.set_tight_layout(True)
	fig1.add_subplot(221)
	
	pl.semilogx(input_sizes, throughputs_MiB, 
				'bo-', ms=MARKER_SIZE, mew=0, mec='b')
	pl.xlabel("fixed file size (Bytes)")
	pl.ylabel("throughput (MiB/sec)")
	pl.text(2E3, 27, "(A)")

	fig1.add_subplot(222)
	pl.loglog(input_sizes, resp_times, 
			  'mo-', ms=MARKER_SIZE, mew=0, mec='m')
	pl.xlabel("fixed file size (Bytes)")
	pl.ylabel("response time (sec)")
	pl.text(2E3, 500, "(B)")

	fig1.add_subplot(223)
	pl.semilogx(resp_times, throughputs_MiB, 
				'go-', ms=MARKER_SIZE, mew=0, mec='g')
	pl.xlabel("response time(sec)")
	pl.ylabel("throughput (MiB/sec)")
	pl.text(0.2, 27, "(C)")

	pl.tight_layout()
	pl.savefig("./figures/%s" % result_file.replace(".txt", ".pdf"))

def bode(G,f=np.arange(.01,100,.01)):
    plt.figure()
    jw = 2*np.pi*f*1j
    y = np.polyval(G.num, jw) / np.polyval(G.den, jw)
    mag = 20.0*np.log10(abs(y))
    phase = np.arctan2(y.imag, y.real)*180.0/np.pi % 360

    plt.subplot(211)
    #plt.semilogx(jw.imag, mag)
    plt.semilogx(f,mag)
    plt.grid()
    plt.gca().xaxis.grid(True, which='minor')

    plt.ylabel(r'Magnitude (db)')

    plt.subplot(212)
    #plt.semilogx(jw.imag, phase)
    plt.semilogx(f,phase)
    plt.grid()
    plt.gca().xaxis.grid(True, which='minor')
    plt.ylabel(r'Phase (deg)')
    plt.yticks(np.arange(0, phase.min()-30, -30))

    return mag, phase

def createAndSaveFig(xData, yData, figFileRoot, xLabel="", yLabel="",
                     fileExt='.png', xMin=-1, xMax=-1, yMin=-10, yMax=-1,
                     log2Y=0, log2X=0, plotType='bo', axisFontSize=20,
                     tickFontSize=16, svgFlag=0):

    figFileName = figFileRoot + fileExt
    xData = convert_list_to_array(xData)
    yData = convert_list_to_array(yData)
    if log2Y == 1:
        tempPlot = plt.semilogy(xData, yData, plotType, basey=2, hold='False')
    else:
        if log2X == 0:
            tempPlot = plt.plot(
                xData, yData, plotType, hold="False", alpha=0.5)
        else:
            tempPlot = plt.semilogx(
                xData, yData, plotType, basey=2, hold='False')
    plt.xlabel(xLabel, fontsize=axisFontSize)
    plt.ylabel(yLabel, fontsize=axisFontSize)
    ax = plt.gca()
    for tick in ax.xaxis.get_major_ticks():
        tick.label1.set_fontsize(tickFontSize)
    for tick in ax.yaxis.get_major_ticks():
        tick.label1.set_fontsize(tickFontSize)
    if xMin == -1:
        xMin = min(xData.tolist())
    if xMax == -1:
        xMax = max(xData.tolist())
    if yMin == -10:
        yMin = min(yData.tolist())
        if isnan(yMin):
            yMin = 0
    if yMax == -1:
        yMax = 0
        yDataList = yData.tolist()
        for item in yDataList:
            if not isnan(item):
                if item > yMax:
                    yMax = item
    plt.xlim(xMin, xMax)
    plt.ylim(yMin, yMax)
    plt.savefig(figFileName, dpi=150)
    if svgFlag == 1:
        figFileName = figFileRoot + '.svg'
        plt.savefig(figFileName, dpi=150)
    plt.clf()

    print device.ask("SSTR ? 0")
    print device.ask("SSTP ? 0")

    device.write("SRPT 2,0")

    device.ask("DSPS?")
    device.write("STRT")
    dataa, datab = False, False
    while not dataa or not datab:
        res = device.display_status_word(int(device.ask("DSPS ?")))
        print res
        codes = res[1]
        if 'SSA' in codes:
            dataa = True
        if 'SSB' in codes:
            datab = True
        # if dataa and datab:
        #     break
        sleep(1)

    # print device.ask("ACTD ?")
    # print device.ask("DTRD ? 2,0")
    # device.write("NOTE 0,1,0,50,50,Hello")
    # print device.ask("DUMP")
    data = [float(num) for num in device.ask("DSPY ? 0").split(',')]
    data = data[0:-1]
    pts = len(data)
    f = np.logspace(np.log10(start), np.log10(stop), pts) # this is incorrect
    pl.semilogx(f, data)
    pl.show()