def plot_xy(self, x, y, xerror=[], yerror=[], title=' ', xLabel=' ', yLabel=' '):
        """
        Simple X vs Y plot

        Inputs:
        ------
          - x = 1D array
          - y = 1D array

        Keywords:
        --------
          - xerror = error on 'x', 1D array
          - yerror = error on 'y', 1D array
          - title = plot title, string
          - xLabel = title of the x-axis, string
          - yLabel = title of the y-axis, string
        """
        fig = plt.figure(figsize=(18,10))
        plt.rc('font',size='22')
        self._fig = plt.plot(x, y, label=title)
        scale = 1
        ticks = ticker.FuncFormatter(lambda lon, pos: '{0:g}'.format(lon/scale))
        plt.ylabel(yLabel)
        plt.xlabel(xLabel)
        if not yerror==[]:
            #plt.errorbar(x, y, yerr=yerror, fmt='o', ecolor='k')
            plt.fill_between(x, y-yerror, y+yerror,
            alpha=0.2, edgecolor='#1B2ACC', facecolor='#089FFF', antialiased=True)
        if not xerror==[]:
            #plt.errorbar(x, y, xerr=xerror, fmt='o', ecolor='k')
            plt.fill_betweenx(y, x-xerror, x+xerror,
            alpha=0.2, edgecolor='#1B2ACC', facecolor='#089FFF', antialiased=True)

        plt.show() 

def plot_wiggle(self,figsize=[5,10],fill=True,perc=100,scale=1,subplot=False):
	# plot wiggle traces
	# figsize=[5,10]: matplotlib figure size [inch]
	# fill=True: fill values greater than zero
	# perc=100: percent clip
	# scale=1: scale trace plots
	plotdata=self.data
	plotdata=perc_clip(plotdata,perc)
	print("min=%s max=%s"%(plotdata.min(),plotdata.max()))
	maxval=np.abs(plotdata).max()
	ns=pk.get_ns(self)
	dt=pk.get_dt(self)
	ntr=pk.get_ntr(self)
	t=np.arange(ns)*dt

	if not subplot: plt.figure(figsize=figsize)
	for itr in range(ntr):
		trace=plotdata[itr,:]
		x=itr+trace/maxval*scale
		plt.plot(x,t,'k-',linewidth=0.5)
		if fill: plt.fill_betweenx(t,x,itr,where=x>itr,color='black',linewidth=0.)

	plt.xlim([-2,ntr+1])
	plt.ylim([t[-1],t[0]])
	plt.gca().xaxis.tick_top()
	plt.gca().xaxis.set_label_position('top')
	if subplot: return
	plt.xlabel('Trace number',fontsize='large')
	plt.ylabel('Time (s)',fontsize='large')

def plot_vawig(axhdl, data, t, excursion):

    import numpy as np
    import matplotlib.pyplot as plt

    [ntrc, nsamp] = data.shape
    

    
    
    t = np.hstack([0, t, t.max()])
    
    for i in range(0, ntrc):
        tbuf = excursion * data[i,:] / np.max(np.abs(data)) + i
        
        tbuf = np.hstack([i, tbuf, i])
            
        axhdl.plot(tbuf, t, color='black', linewidth=0.5)
        plt.fill_betweenx(t, tbuf, i, where=tbuf>i, facecolor=[0.6,0.6,1.0], linewidth=0)
        plt.fill_betweenx(t, tbuf, i, where=tbuf<i, facecolor=[1.0,0.7,0.7], linewidth=0)
    
    axhdl.set_xlim((-excursion, ntrc+excursion))
    axhdl.xaxis.tick_top()
    axhdl.xaxis.set_label_position('top')
    axhdl.invert_yaxis()

def plot_fill(llr_cur, tkey, asimov_llr, hist_vals, bincen, fit_gauss, **kwargs):
    """
    Plots fill between the asimov llr value and the histogram values
    which represent an LLR distribution.
    """
    validate_key(tkey)

    expr = 'bincen < asimov_llr' if 'true_N' in tkey else 'bincen > asimov_llr'

    plt.fill_betweenx(
        hist_vals, bincen, x2=asimov_llr, where=eval(expr), **kwargs)

    pvalue = (1.0 - float(np.sum(llr_cur > asimov_llr))/len(llr_cur)
              if 'true_N' in tkey else
              (1.0 - float(np.sum(llr_cur < asimov_llr))/len(llr_cur)))

    sigma_fit = np.fabs(asimov_llr - fit_gauss[1])/fit_gauss[2]
    #logging.info(
    #    "  For tkey: %s, gaussian computed mean (of alt MH): %.3f and sigma: %.3f"
    #    %(tkey,fit_gauss[1],fit_gauss[2]))
    pval_gauss = 1.0 - norm.cdf(sigma_fit)
    sigma_1side = np.sqrt(2.0)*erfinv(1.0 - pval_gauss)

    mctrue_row = [tkey,asimov_llr,llr_cur.mean(),pvalue,pval_gauss,sigma_fit,
                  sigma_1side]

    return mctrue_row

def plot_posteriors(model_params,plotfile,nbins=30,names=None):

    import matplotlib
    import matplotlib.pyplot as plt

    #----get array dimensions (number of parameters)----
    npar = model_params.shape[1]

    #----loop through parameters----
    for p in xrange(npar):
        
        #--skip iteration column or if parameter is fixed--
        if (max(model_params[:,p]) != min(model_params[:,p])) and p != 0:
            
            y = model_params[:,p]        
            y_hist, x_bin = np.histogram(y,bins=nbins)
            fig,ax=plt.subplots()
            plt.bar(x_bin[:-1],y_hist,width=x_bin[1]-x_bin[0])
            if names != None:
                plt.xlabel(names[p])
            ymin,ymax = ax.get_ylim()

            #-plot median and 1sigma range-
            med = np.median(y)
            sig = np.std(y)
            plt.plot([med,med],[ymin,ymax],color='red',linestyle='-')
            plt.fill_betweenx([0.0,ymax],[med-sig,med-sig],[med+sig,med+sig],color='red',alpha=0.2)

            #-save plot-
            plotfile.savefig()
            plt.close()

def make_group_cumdist_fig(LD,group_by,pcol):
    '''Make a plot showing the cumulative distribution of the within-group avg of a 
    specified response variable. compared to the cum dist expected by chance.
    INPUTS: 
        LD: pandas dataframe
        group_by: column to groupby for computing avgs.
        pcol: name of response variable column.
    OUTPUTS:
        fig: figure handle'''    

    name_legend_map = {'counts': 'Number of loans (thousands)',
				 'ROI': 'Average ROI (%)',
				  'int_rate': 'interest rate (%)',
				  'default_prob': 'default probability',
				  'dti': 'Debt-to-income ratio',
				  'emp_length': 'employment length (months)',
                        'annual_inc': 'annual income ($)'}

    min_group_loans = 100 #only use states with at least this many loans
    good_groups = LD.groupby(group_by).filter(lambda x: x[pcol].count() >= min_group_loans)
    n_groups = len(good_groups[group_by].unique())
    group_stats = good_groups.groupby(group_by)[pcol].agg(['mean','sem'])
    group_stats.sort_values(by='mean',inplace=True)
    ov_avg = good_groups[pcol].mean()
    
    #compute bootstrap estimates of null distribution of group-avgs
    boot_samps = 500 #number of bootstrap samples to use when estimating null dist
    shuff_avgs = np.zeros((boot_samps,n_groups))
    shuff_data = good_groups.copy()
    for cnt in xrange(boot_samps):
        shuff_data[pcol] = np.random.permutation(shuff_data[pcol].values)
        shuff_avgs[cnt,:] = shuff_data.groupby(group_by)[pcol].mean().sort_values()
    
    yax = np.arange(n_groups)
    rmean = np.mean(shuff_avgs,axis=0)
    rsem = np.std(shuff_avgs,axis=0)

    #plot avg and SEM of within-state returns
    fig, ax1 = plt.subplots(1,1,figsize=(6.0,5.0))
    if group_by == 'zip3':
        ax1.errorbar(group_stats['mean'].values,yax)
    else:
        ax1.errorbar(group_stats['mean'],yax,xerr=group_stats['sem'])
        plt.fill_betweenx(yax, rmean-rsem, rmean+rsem,facecolor='r',alpha=0.5,linewidth=0)
        plt.yticks(yax, group_stats.index,fontsize=6)
    
    ax1.plot(rmean,yax,'r')
    plt.legend(['Measured-avgs','Shuffled-avgs'],loc='best')
    if group_by == 'zip3':    
        plt.ylabel('Zip codes',fontsize=16)
    else:    
        plt.ylabel('States',fontsize=16)
    plt.xlabel(name_legend_map[pcol],fontsize=16)   
    plt.xlim(np.min(group_stats['mean']),np.max(group_stats['mean']))
    plt.ylim(0,n_groups)
    ax1.axvline(ov_avg,color='k',ls='dashed')
    plt.tight_layout()

    return fig

def plot_mutation_rate_violins(libraries, out_prefix, nucleotides_to_count='ATCG', exclude_constitutive=False):
    #Makes violin plots of raw mutation rates
    data = []
    labels = []
    for library in libraries:
        labels.append(library.lib_settings.sample_name)
        data.append([math.log10(val) for val in library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                        nucleotides_to_count=nucleotides_to_count,
                                                        exclude_constitutive=exclude_constitutive) if val>0])

    colormap = uniform_colormaps.viridis
    fig = plt.figure(figsize=(5,8))
    ax1 = fig.add_subplot(111)

    # Hide the grid behind plot objects
    ax1.yaxis.grid(True, linestyle='-', which='major', color='lightgrey', alpha=0.5)
    ax1.set_axisbelow(True)

    #ax1.set_xlabel(ylabel)
    plt.subplots_adjust(left=0.1, right=0.95, top=0.9, bottom=0.25)

    pos = range(1,len(libraries)+1)  # starts at 1 to play nice with boxplot
    dist = max(pos)-min(pos)
    w = min(0.15*max(dist,1.0),0.5)
    for library,p in zip(libraries,pos):
        d = [math.log10(val) for val in library.list_mutation_rates(subtract_background=False, subtract_control=False,
                                                        nucleotides_to_count=nucleotides_to_count,
                                                        exclude_constitutive=exclude_constitutive) if val>0]
        k = stats.gaussian_kde(d) #calculates the kernel density
        m = k.dataset.min() #lower bound of violin
        M = k.dataset.max() #upper bound of violin
        x = numpy.arange(m,M,(M-m)/100.) # support for violin
        v = k.evaluate(x) #violin profile (density curve)
        v = v/v.max()*w #scaling the violin to the available space
        plt.fill_betweenx(x,p,v+p,facecolor=colormap((p-1)/float(len(libraries))),alpha=0.3)
        plt.fill_betweenx(x,p,-v+p,facecolor=colormap((p-1)/float(len(libraries))),alpha=0.3)
    if True:
        bplot = plt.boxplot(data,notch=1)
        plt.setp(bplot['boxes'], color='black')
        plt.setp(bplot['whiskers'], color='black')
        plt.setp(bplot['fliers'], color='red', marker='.')

    per50s = []
    i = 1
    for datum in data:
        #per50s.append(stats.scoreatpercentile(datum, 50))
        t = stats.scoreatpercentile(datum, 50)

        per50s.append(t)
        #ax1.annotate(str(round(t,3)), xy=(i+0.1, t), xycoords='data', arrowprops=None, fontsize='small', color='black')
        i+= 1
    #ax1.set_xticks([0.0, 0.5, 1.0, 1.5])
    #ax1.set_yscale('log')
    ax1.set_ylabel('log10 mutation rate')
    ax1.set_ylim(-5, 0)
    xtickNames = plt.setp(ax1, xticklabels=labels)
    plt.setp(xtickNames, rotation=90, fontsize=6)
    plt.savefig(out_prefix+'_logviolin.pdf', transparent='True', format='pdf')
    plt.clf()

def seismic_wiggle(section, dt=0.004, ranges=None, scale=1.,
                   color='k', normalize=False):
    """
    Plot a seismic section (numpy 2D array matrix) as wiggles.

    Parameters:

    * section :  2D array
        matrix of traces (first dimension time, second dimension traces)
    * dt : float
        sample rate in seconds (default 4 ms)
    * ranges : (x1, x2)
        min and max horizontal values (default trace number)
    * scale : float
        scale factor multiplied by the section values before plotting
    * color : tuple of strings
        Color for filling the wiggle, positive  and negative lobes.
    * normalize :
        True to normalizes all trace in the section using global max/min
        data will be in the range (-0.5, 0.5) zero centered

    .. warning::
        Slow for more than 200 traces, in this case decimate your
        data or use ``seismic_image``.

    """
    npts, ntraces = section.shape  # time/traces
    if ntraces < 1:
        raise IndexError("Nothing to plot")
    if npts < 1:
        raise IndexError("Nothing to plot")
    t = numpy.linspace(0, dt*npts, npts)
    amp = 1.  # normalization factor
    gmin = 0.  # global minimum
    toffset = 0.  # offset in time to make 0 centered
    if normalize:
        gmax = section.max()
        gmin = section.min()
        amp = (gmax-gmin)
        toffset = 0.5
    pyplot.ylim(max(t), 0)
    if ranges is None:
        ranges = (0, ntraces)
    x0, x1 = ranges
    # horizontal increment
    dx = float((x1-x0)/ntraces)
    pyplot.xlim(x0, x1)
    for i, trace in enumerate(section.transpose()):
        tr = (((trace-gmin)/amp)-toffset)*scale*dx
        x = x0+i*dx  # x positon for this trace
        pyplot.plot(x+tr, t, 'k')
        pyplot.fill_betweenx(t, x+tr, x, tr > 0, color=color)

def drawhigh(cid,filesize,view,threshold,high,lowlimit=0,highlimit=0):
    avgview=sum(view[5:-5])/len(view)
    highdur=map(lambda x:(x[0],x[-1]),high)

    # 图像设置
    plt.figure(figsize=(15,7)) # figsize()设置的宽高比例是是15:7，图片的尺寸会根据这个比例进行调节
    #plt.xlim(-3,19)
    ylow=min(view)-500      #y轴下限
    yhigh=max(view)+500     #y轴上限
    plt.ylim(ylow,yhigh)
    plt.grid(which='both')


    #绘制结果数据
    plt.plot(range(1,len(view)+1),view,'bo-',ms=1,lw=0.5,label='origin')      # 原始图像
    if lowlimit and highlimit:
        plt.axhline(y=lowlimit,lw=3,ls='-',color='m',label='lowlimit')        # 低限
        plt.axhline(y=highlimit,lw=3,ls='-',color='m',label='highlimit')      # 高限

    plt.axhline(y=avgview,lw=1,ls='--',color='b',label='mean')                # 均值
    #plt.axhline(y=avgview*adjustv,lw=1,ls='--',color='g',label='mean*1.2')   # 均值*adjustv
    plt.axhline(y=threshold,lw=2,ls='--',color='r',label='threshold=1.2*mean')       # 阈值
    plt.legend(loc='upper right')
    plt.xlabel('time (s)')
    plt.ylabel('views')

    # 标注高潮区间
    for item in highdur:
        #plt.axvline(x=item[0],lw=2)
        #plt.axvline(x=item[1],lw=2)
        plt.annotate('',xy=(item[1],threshold),xytext=(item[0],threshold),arrowprops=dict(arrowstyle="->",connectionstyle="arc3",color='g'))
        plt.fill_betweenx([ylow,yhigh],item[0], item[1], linewidth=1, alpha=0.2, color='r')

    plt.show()

    # 结果保存
    '''
    resultpath='D:\\hot_pic2'
    if not os.path.exists(resultpath):
        os.mkdir(resultpath)
    fname=os.path.join(resultpath,cid+'.'+str(filesize)+'.jpg')
    print fname
    plt.savefig(fname,dpi = 300)
    plt.close()
    '''
    return 0;

    def show_distr(x,y, vertical=False, label=None, color='blue', linecolor='k', quantile=False):
        var_2d=np.copy(y)
        if vertical:
            var_2d=np.copy(x)
        mid=np.nanmean(var_2d, axis=0)
        lower=mid - np.nanstd(var_2d, axis=0)
        upper=mid + np.nanstd(var_2d, axis=0)
        if quantile:
            lower=np.nanpercentile(var_2d,25,axis=0)
            upper=np.nanpercentile(var_2d,75,axis=0)

        if vertical:
            plt.fill_betweenx(y,lower,upper, color=color)
            plt.plot(mid,y,color=linecolor, linewidth=2)
        else:
            plt.fill_between(x,lower,upper, color=color)
            plt.plot(x,mid,color=linecolor, linewidth=2)

def plot_airmass(objectlist,obsvat,date):
    #ax = plt.subplot(111)
    
    for obj in objectlist:
        altdata = compute_alt_plot(obj[0],obj[1],obsvat,date)
        plt.plot(altdata[:,0],altdata[:,1])


    morn_twilight,even_twilight = functions.calc_twilight(obsvat,date)
    

    plt.fill_betweenx([0,90],[morn_twilight.datetime(),morn_twilight.datetime()],x2=[even_twilight.datetime(),even_twilight.datetime()],color="0.5")

    locs,labels = plt.xticks()
    plt.setp(labels,rotation=45)
    plt.ylim(0,90)
    plt.show()

 def peek(self):
     plt.imshow(self.data, aspect='auto', interpolation='none',
                extent=[self.times[0].to(u.s).value,
                        self.times[-1].to(u.s).value,
                        self.latitude[0].to(u.degree).value,
                        self.latitude[-1].to(u.degree).value])
     plt.xlim(0, self.times[-1].to(u.s).value)
     if self.times[0].to(u.s).value > 0.0:
         plt.fill_betweenx([self.latitude[0].to(u.degree).value,
                            self.latitude[-1].to(u.degree).value],
                           self.times[0].to(u.s).value,
                           hatch='X', facecolor='w', label='not observed')
     plt.ylabel('degrees of arc from first measurement')
     plt.xlabel('time since originating event (seconds)')
     plt.title('arc: ' + self.title)
     plt.legend(framealpha=0.5)
     plt.show()
     return None

def plotBootROC(rocDfL, labelL=None, aucL=None, ciParam='fpr'):
    """Plot of ROC curves with confidence intervals.

    Parameters
    ----------
    rocDfL : list of pd.DataFrames
        Each DataFram is one model and must include columns
        fpr_est, tpr_est, fpr_lb, fpr_ub
    labelL : list of str
        Names of each model for legend
    aucL : list of floats
        AUC scores of each model for legend"""
    if labelL is None and aucL is None:
        labelL = ['Model %d' % i for i in range(len(rocDfL))]
    elif labelL is None:
        labelL = ['Model %d (AUC = %0.2f [%0.2f, %0.2f])' % (i, auc[0], auc[1], auc[2]) for i, auc in enumerate(aucL)]
    else:
        labelL = ['%s (AUC = %0.2f [%0.2f, %0.2f])' % (label, auc[0], auc[1], auc[2]) for label, auc in zip(labelL, aucL)]

    colors = sns.color_palette('Set1', n_colors=len(rocDfL))

    plt.cla()
    plt.gca().set_aspect('equal')
    for i, (rocDf, label) in enumerate(zip(rocDfL, labelL)):
        if ciParam == 'fpr':
            plt.fill_betweenx(rocDf['tpr_est'], rocDf['fpr_lb'], rocDf['fpr_ub'], alpha=0.3, color=colors[i])
        elif ciParam == 'tpr':
            plt.fill_between(rocDf['fpr_est'], rocDf['tpr_lb'], rocDf['tpr_ub'], alpha=0.3, color=colors[i])
        plt.plot(rocDf['fpr_est'], rocDf['tpr_est'],'-', color=colors[i], lw=2)
        # plt.plot(rocDf['fpr_est'], rocDf['tpr_lb'], '.--', color=colors[i], lw=1)
        # plt.plot(rocDf['fpr_est'], rocDf['tpr_ub'], '.--', color=colors[i], lw=1)
        # plt.plot(rocDf['fpr_lb'], rocDf['tpr_est'], '--', color=colors[i], lw=1)
        # plt.plot(rocDf['fpr_ub'], rocDf['tpr_est'], '--', color=colors[i], lw=1)
    plt.plot([0, 1], [0, 1], '--', color='gray', label='Chance')
    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('ROC')
    plt.legend([plt.Line2D([0, 1], [0, 1], color=c, lw=2) for c in colors], labelL, loc='lower right', fontsize=10)
    plt.show()

def gen_fill(tail_num, deg_left, deg_right):
    if tail_num == 3:
        plt.fill_betweenx(mlab.normpdf(x,mean,sigma),deg_left,x, where = ( x <= deg_left))
        plt.fill_betweenx(mlab.normpdf(x,mean,sigma), x,deg_right, where = ( x >= deg_right))    
        plt.draw()
        plt.show()
    elif tail_num ==  2: #right
        plt.fill_betweenx(mlab.normpdf(x,mean,sigma), x,deg_right, where = ( x >= deg_right))    
        plt.draw()
        plt.show()
    elif tail_num == 1 :
        plt.fill_betweenx(mlab.normpdf(x,mean,sigma),deg_left,x, where = ( x <= deg_left))
        plt.draw()
        plt.show()

def plot_wiggle(data,zz=1,skip=1,gain=1,alpha=0.7,black=False):
    '''
    Wiggle plot of generic 2D numpy array.

    INPUT
    data: 2D numpy array
    zz: vertical sample rate in depth or time
    skip: interval to choose traces to draw
    gain: multiplier applied to each trace
    '''
    [n_samples,n_traces]=data.shape
    t=range(n_samples)
    plt.figure(figsize=(9.6,6))
    for i in range(0, n_traces,skip):
        trace=gain*data[:,i] / np.max(np.abs(data))
        plt.plot(i+trace,t,color='k', linewidth=0.5)
        if black==False:
            plt.fill_betweenx(t,trace+i,i, where=trace+i>i, facecolor=[0.6,0.6,1.0], linewidth=0)
            plt.fill_betweenx(t,trace+i,i, where=trace+i<i, facecolor=[1.0,0.7,0.7], linewidth=0)
        else:
            plt.fill_betweenx(t,trace+i,i, where=trace+i>i, facecolor='black', linewidth=0, alpha=alpha)
    locs,labels=plt.yticks()
    plt.yticks(locs,[n*zz for n in locs.tolist()])
    plt.grid()
    plt.gca().invert_yaxis()

def plot_dDAR():
	lam_mu = 0.3+np.arange(2200)/1000
	lam_mu_ref = 0.55
	zeta_deg = 10*np.arange(8)
	plt.ylim([-3,5])
	##
	## mark X-SHOOTER arms
	plt.fill_betweenx(plt.ylim(),0.3,0.55,color="blue",alpha=0.3)
	plt.fill_betweenx(plt.ylim(),0.55,1.0,color="green",alpha=0.3)
	plt.fill_betweenx(plt.ylim(),1.0,2.5,color="red",alpha=0.3)

	for z in zeta_deg:
		am = 1/np.sin(np.deg2rad(90-z))
		plt.plot(lam_mu, dDAR(z, lam_mu, lam_mu_ref), label="{0} deg, airmass {1:5.2f}".format(z,am))
	
	plt.legend()
	plt.title("Differential atmospheric refraction w.r.t. {0:5.2f} micron".format(lam_mu_ref))
	plt.xlabel("Wavelength in micron")
	plt.ylabel("Differential atmospheric refraction in arcsec")

	plt.text(2.4, -2.5, "Paranal standard conditions: 11.5 deg C, 743 hPa, 14.5 % RH", fontsize=9, ha="right")

import matplotlib.pyplot as plt
inds_t, inds_v = [(clusters[cluster_ind]) for ii, cluster_ind in
                  enumerate(good_cluster_inds)][0]  # first cluster

times = np.arange(X[0].shape[1]) * tstep * 1e3

plt.clf()
colors = ['y', 'b', 'g', 'purple']
event_ids = ['l_aud', 'r_aud', 'l_vis', 'r_vis']

for ii, (condition, color, eve_id) in enumerate(zip(X, colors, event_ids)):
    # extract time course at cluster vertices
    condition = condition[:, :, inds_v]
    # normally we would normalize values across subjects but
    # here we use data from the same subject so we're good to just
    # create average time series across subjects and vertices.
    mean_tc = condition.mean(axis=2).mean(axis=0)
    std_tc = condition.std(axis=2).std(axis=0)
    plt.plot(times, mean_tc.T, color=color, label=eve_id)
    plt.fill_between(times, mean_tc + std_tc, mean_tc - std_tc, color='gray',
                     alpha=0.5, label='')

plt.xlabel('Time (ms)')
plt.ylabel('Activation (F-values)')
plt.xlim(times[[0, -1]])
plt.fill_betweenx(np.arange(*plt.ylim()), times[inds_t[0]],
                  times[inds_t[-1]], color='orange', alpha=0.3)
plt.legend()
plt.title('Interaction between stimulus-modality and location.')
plt.show()

 def horiz_plot(v_coord, orography, style_args):
     y = v_coord.points
     x = orography.points
     return plt.fill_betweenx(y, x, **style_args)

cb.set_label('Brightness contrast, log10(Shell / BG)')
plt.xlabel('Projected distance from Trapezium, D / arcmin')
plt.ylabel('Bowshock radius, r0 / arcsec')
plt.text(0.05, 0.05, 'Symbol size indicates shell relative thickness, H',
         transform=ax.transAxes, fontsize='x-small')
ax.set_xlim(0.05, 20.0)
ax.set_ylim(0.3, 11.0)
ax.set_xscale('log')
ax.set_yscale('log')
fig.savefig(pltfile)
figlist.append('[[file:luis-programas/{0}][{0}]]'.format(pltfile))

pltfile = 'will-PA-vs-PA.pdf'
fig = plt.figure(figsize=(7,6))
ax = fig.add_subplot(111, axisbg="#eeeeee")
plt.fill_betweenx([-90.0, 90.0], [0.0, 0.0], [90.0, 90.0], zorder=-10, alpha=0.05)
plt.fill_betweenx([-90.0, 90.0], [180.0, 180.0], [270.0, 270.0], zorder=-10, alpha=0.05)
plt.fill_betweenx([-90.0, 90.0], [360.0, 360.0], [450.0, 450.0], zorder=-10, alpha=0.05)
plt.text(45.0, -80.0, 'NE\nquadrant',  ha='center', fontsize='x-small')
plt.text(135.0, -80.0, 'SE\nquadrant', ha='center', fontsize='x-small')
plt.text(225.0, -80.0, 'SW\nquadrant', ha='center', fontsize='x-small')
plt.text(315.0, -80.0, 'NW\nquadrant', ha='center', fontsize='x-small')
plt.axhline(zorder=-5)
plt.scatter(PA_star[m], dPA[m], s=20*tab['R_out'][m], c=D60[m], cmap=plt.cm.hot, alpha=0.6)
label_sources(tab['Object'], PA_star, dPA, np.abs(dPA) > 45.0, allmask=m)
cb = plt.colorbar()
cb.set_label('Projected distance from Trapezium, D / arcmin')
plt.xlabel('PA of source from Trapezium, deg')
plt.ylabel('Angle between bowshock axis and radial direction, deg')
ax.set_xlim(-30.0, 375.0)
ax.set_ylim(-90.0, 90.0)

def plot_llr_distributions(llr_nmh, llr_imh, nbins, plot_gauss=True, imh_true=False):
    """Plots LLR distributions-expects llr_nmh and llr_imh to be type
    Series. Also plots vertical line signifying the mean of the
    hierarchy assumed to be given in nature, and the percentage of
    trials from the opposite hierarchy with LLR beyond the mean.
    """

    fig = plt.figure(figsize=(9, 8))
    label_text = r"$\log ( \mathcal{L}(data: IH|IH) / \mathcal{L}( data: IH|NH) )$"
    llr_imh.hist(bins=nbins, histtype="step", lw=2, color="b", label=label_text)
    hist_vals_imh, bincen_imh = plot_error(llr_imh, nbins, fmt=".b", lw=2)
    if plot_gauss:
        fit_imh = plot_gauss_fit(llr_imh, hist_vals_imh, bincen_imh, color="b", lw=2)

    label_text = r"$\log ( \mathcal{L}(data: NH|IH) / \mathcal{L}(data: NH|NH) )$"
    llr_nmh.hist(bins=nbins, histtype="step", lw=2, color="r", label=label_text)
    hist_vals_nmh, bincen_nmh = plot_error(llr_nmh, nbins, fmt=".r", lw=2)
    if plot_gauss:
        fit_nmh = plot_gauss_fit(llr_nmh, hist_vals_nmh, bincen_nmh, color="r", lw=2)

    if imh_true:
        mean_val = llr_nmh.mean()
        pvalue = 1.0 - float(np.sum(llr_imh > mean_val)) / len(llr_imh)
    else:
        mean_val = llr_imh.mean()
        pvalue = float(np.sum(llr_nmh > mean_val)) / len(llr_nmh)

    ymax = max(hist_vals_nmh) if imh_true else max(hist_vals_imh)
    bincen = bincen_imh if imh_true else bincen_nmh
    vline = plt.vlines(mean_val, 1, ymax, colors="k", linewidth=2, label=("pval = %.4f" % pvalue))

    sigma_1side = np.sqrt(2.0) * erfinv(1.0 - pvalue)
    sigma_2side = norm.isf(pvalue)
    print "  Using non-gauss fit: "
    print "    pvalue: %.5f" % pvalue
    print "    sigma 1 sided (erfinv): %.4f" % sigma_1side
    print "    sigma 2 sided (isf)   : %.4f" % sigma_2side

    sigma_fit_imh = (fit_imh[1] - fit_nmh[1]) / fit_imh[2]
    sigma_fit_nmh = (fit_imh[1] - fit_nmh[1]) / fit_nmh[2]
    pval_imh = 1.0 - norm.cdf(sigma_fit_imh)
    pval_nmh = 1.0 - norm.cdf(sigma_fit_nmh)
    sigma_1side_imh = np.sqrt(2.0) * erfinv(1.0 - pval_imh)
    sigma_1side_nmh = np.sqrt(2.0) * erfinv(1.0 - pval_nmh)

    print "\n\n  pval IMH: %.5f, pval NMH: %.5f" % (pval_imh, pval_nmh)
    print "  sigma gauss fit (IMH true): %.4f/%.4f" % (sigma_fit_imh, sigma_1side_imh)
    print "  sigma gauss fit (NMH true): %.4f/%.4f" % (sigma_fit_nmh, sigma_1side_nmh)

    sigma_error_nmh = np.sqrt(1.0 + 0.5 * (2.0 / sigma_1side_nmh) ** 2) / np.sqrt(len(llr_nmh))
    sigma_error_imh = np.sqrt(1.0 + 0.5 * (2.0 / sigma_1side_imh) ** 2) / np.sqrt(len(llr_imh))
    logging.info("total trials: %d", len(llr_nmh))
    logging.info("  nmh sigma error: %f" % sigma_error_nmh)
    logging.info("  imh sigma error: %f" % sigma_error_imh)

    if imh_true:
        plt.fill_betweenx(hist_vals_imh, bincen, x2=mean_val, where=bincen < mean_val, alpha=0.5, hatch="xx")
    else:
        plt.fill_betweenx(hist_vals_nmh, bincen, x2=mean_val, where=bincen > mean_val, alpha=0.5, hatch="xx")

    if args.present:
        plt.ylabel("# Trials")
        plt.xlabel("LLR value")
    else:
        plt.ylabel("# Trials", fontsize="x-large")
        plt.xlabel("LLR value", fontsize="x-large")

    return