def plot_spectrograms(bsl,rec,rate,title):
    plt.close()
    fig, ax = plt.subplots(nrows=9, ncols=2, sharex='col', sharey='row')
    plt.subplots_adjust(wspace = .05,hspace = 0.4 )
    ny_nfft=1024
    i=0
    while i<9:
        Pxx, freq, bins, im = ax[i,0].specgram(bsl[i],NFFT=ny_nfft,Fs=rate)
        ax[i,0].set_ylim([0, 40])
        if(i==8):
            ax[i,0].set_xlabel("Time, seconds")
        ax[i,0].set_ylabel("Freq, Hz")
        ax[i,0].set_title(title+' baseline sleep, REM stage, Channel:'+str(i+1))
        i=i+1
    i=0
    while i<9:
        Pxx, freq, bins, im = ax[i,1].specgram(rec[i],NFFT=ny_nfft,Fs=rate)
        #ax[i,1].ylim(0,40)
        ax[i,1].set_ylim([0, 40])
        #ax[i,1].set_xlim([0, 10000]) #13000])
        if(i==8):
            ax[i,1].set_xlabel("Time, seconds")
        #ax[i,1].set_ylabel("Freq, Hz")
        ax[i,1].set_title(title+'   recovery sleep, REM stage, Channel:'+str(i+1))
        i=i+1
    plt.show()
    return

def plot_setup_post(figure_number = None, show = True, save_file = None,
					legend = True, legend_location = 0):
	"""
	Handles post-figure setup, including legends, file saving (save_file is
	desired filename), showing the figure, and clearing it.
	"""

	if figure_number:
		pyp.figure(figure_number)

	pyp.subplots_adjust(bottom=.5) # adjustment to give more xlabel space

	# change limits
	# pyp.xlim( xmin = 0, xmax = 10000 )
	# pyp.ylim( ymin = 0, ymax = 10000 )
	# pyp.ylim( (0,10000) ) # equivalent to the line above

	if legend:
		pyp.legend(loc = legend_location)
	if save_file:
		pyp.savefig(save_file)
	if show:
		pyp.show()
	else:
		pyp.clf() # clears figure if not plotted

def plot_data(tag):
    data_array = tag.references[0]
    voltage = np.zeros(data_array.data.shape)
    data_array.data.read_direct(voltage)

    x_axis = data_array.dimensions[0]
    time = x_axis.axis(data_array.data_extent[0])

    spike_times = tag.positions[:]

    feature_data_array = tag.features[0].data
    snippets = tag.features[0].data[:]
    single_snippet = tag.retrieve_feature_data(3, 0)[:]

    snippet_time_dim = feature_data_array.dimensions[1]
    snippet_time = snippet_time_dim.axis(feature_data_array.data_extent[1])

    response_axis = plt.subplot2grid((2, 2), (0, 0), rowspan=1, colspan=2)
    single_snippet_axis = plt.subplot2grid((2, 2), (1, 0), rowspan=1, colspan=1)
    average_snippet_axis = plt.subplot2grid((2, 2), (1, 1), rowspan=1, colspan=1)

    response_axis.plot(time, voltage, color="dodgerblue", label=data_array.name)
    response_axis.scatter(spike_times, np.ones(spike_times.shape) * np.max(voltage), color="red", label=tag.name)
    response_axis.set_xlabel(x_axis.label + ((" [" + x_axis.unit + "]") if x_axis.unit else ""))
    response_axis.set_ylabel(data_array.label + ((" [" + data_array.unit + "]") if data_array.unit else ""))
    response_axis.set_title(data_array.name)
    response_axis.set_xlim(0, np.max(time))
    response_axis.set_ylim((1.2 * np.min(voltage), 1.2 * np.max(voltage)))
    response_axis.legend()

    single_snippet_axis.plot(snippet_time, single_snippet.T, color="red", label=("snippet No 4"))
    single_snippet_axis.set_xlabel(
        snippet_time_dim.label + ((" [" + snippet_time_dim.unit + "]") if snippet_time_dim.unit else "")
    )
    single_snippet_axis.set_ylabel(
        feature_data_array.label + ((" [" + feature_data_array.unit + "]") if feature_data_array.unit else "")
    )
    single_snippet_axis.set_title("single stimulus snippet")
    single_snippet_axis.set_xlim(np.min(snippet_time), np.max(snippet_time))
    single_snippet_axis.set_ylim((1.2 * np.min(snippets[3, :]), 1.2 * np.max(snippets[3, :])))
    single_snippet_axis.legend()

    mean_snippet = np.mean(snippets, axis=0)
    std_snippet = np.std(snippets, axis=0)
    average_snippet_axis.fill_between(
        snippet_time, mean_snippet + std_snippet, mean_snippet - std_snippet, color="red", alpha=0.5
    )
    average_snippet_axis.plot(snippet_time, mean_snippet, color="red", label=(feature_data_array.name + str(4)))
    average_snippet_axis.set_xlabel(
        snippet_time_dim.label + ((" [" + snippet_time_dim.unit + "]") if snippet_time_dim.unit else "")
    )
    average_snippet_axis.set_ylabel(
        feature_data_array.label + ((" [" + feature_data_array.unit + "]") if feature_data_array.unit else "")
    )
    average_snippet_axis.set_title("spike-triggered average")
    average_snippet_axis.set_xlim(np.min(snippet_time), np.max(snippet_time))
    average_snippet_axis.set_ylim((1.2 * np.min(mean_snippet - std_snippet), 1.2 * np.max(mean_snippet + std_snippet)))

    plt.subplots_adjust(left=0.15, top=0.875, bottom=0.1, right=0.98, hspace=0.35, wspace=0.25)
    plt.show()

def doit():
    # test it out

    L = 1

    R = 10

    theta = np.radians(np.arange(0,361))

    # draw a circle
    plt.plot(R*np.cos(theta), R*np.sin(theta), c="b")


    # draw some people
    angles = [30, 60, 90, 120, 180, 270, 300]

    for l in angles:
        center = ( (R + 0.5*L)*np.cos(np.radians(l)),
                   (R + 0.5*L)*np.sin(np.radians(l)) )
        draw_person(center, L, np.radians(l - 90), color="r")
        L = 1.1*L

    plt.axis("off")

    ax = plt.gca()
    ax.set_aspect("equal", "datalim")


    plt.subplots_adjust(left=0.05, right=0.98, bottom=0.05, top=0.98)
    plt.axis([-1.2*R, 1.2*R, -1.2*R, 1.2*R])

    f = plt.gcf()
    f.set_size_inches(6.0, 6.0)

    plt.savefig("test.png")

def plot_time_domain_waveform(fig, waveform, imag=False, mag=False, 
                              xlim=None, xlabel=r'$tc^3/GM$',
                              ylabel_pol=r'$h_+ + i h_\times$', 
                              ylabel_amp=r'$A$', ylabel_phase=r'$\Phi$', 
                              pol_legend=True, wave_legend=False):
    """Plot the amplitude, phase, and polarizations of a waveform.
    """
    # Polarization plot
    axes = fig.add_subplot(311)
    t = waveform.time
    hcomp = waveform.get_complex()
    label = r'$h_+$' if pol_legend else ''
    line_list = axes.plot(t, hcomp.real, ls='-', label=label)
    color = line_list[0].get_color()
    if imag:
        label = r'$h_\times$' if pol_legend else ''
        axes.plot(t, hcomp.imag, ls='--', c=color, label=label)
    if mag:
        label = r'$|h_+ + ih_\times|$' if pol_legend else ''
        axes.plot(waveform.time, waveform.amp, ls=':', c=color, label=label)
    
    if xlim is not None: axes.set_xlim(xlim)
    axes.set_ylabel(ylabel_pol, fontsize=16)
    axes.set_xticklabels(axes.get_xticks(), fontsize=14)
    axes.set_yticklabels(axes.get_yticks(), fontsize=14)
    axes.minorticks_on()
    axes.tick_params(which='major', width=2, length=8)
    axes.tick_params(which='minor', width=2, length=4)
    axes.xaxis.set_major_formatter(NullFormatter()) # get rid of x-axis numbers
    axes.legend(fontsize=14, loc='best', ncol=3)    
    
    # Amplitude plot
    axes = fig.add_subplot(312)
    axes.plot(waveform.time, waveform.amp, c=color)
    
    if xlim is not None: axes.set_xlim(xlim)
    axes.set_ylabel(ylabel_amp, fontsize=16)
    axes.set_xticklabels(axes.get_xticks(), fontsize=14)
    axes.set_yticklabels(axes.get_yticks(), fontsize=14)
    axes.minorticks_on()
    axes.tick_params(which='major', width=2, length=8)
    axes.tick_params(which='minor', width=2, length=4)
    axes.xaxis.set_major_formatter(NullFormatter()) # get rid of x-axis numbers
    
    # Phase plot
    axes = fig.add_subplot(313)
    label = wave_legend if wave_legend is not False else ''
    axes.plot(waveform.time, waveform.phase, c=color, label=label)

    if xlim is not None: axes.set_xlim(xlim)
    axes.set_xlabel(xlabel, fontsize=16)
    axes.set_ylabel(ylabel_phase, fontsize=16)
    axes.set_xticklabels(axes.get_xticks(), fontsize=14)
    axes.set_yticklabels(axes.get_yticks(), fontsize=14)
    axes.minorticks_on()
    axes.tick_params(which='major', width=2, length=8)
    axes.tick_params(which='minor', width=2, length=4)
    axes.legend(fontsize=14, loc='best', ncol=2)
    
    subplots_adjust(hspace=0.07)

def plot_fits(direction_rates,fit_curve,title):
    """
    This function takes the x-values and the y-values  in units of spikes/s 
    (found in the two columns of direction_rates and fit_curve) and plots the 
    actual values with circles, and the curves as lines in both linear and 
    polar plots.
    """
    #print direction_rates
    
    
    
   # print fit_curve
   
    plt.subplots_adjust(hspace = 0.6)   
    y_max = np.max(direction_rates[:,1]) + 5
    plt.subplot(2,2,3)
    plt.axis([0,360,0,y_max])
    plt.plot(direction_rates[:,0], direction_rates[:,1],'o')
    plt.plot(fit_curve[:,0],fit_curve[:,1], '-')
    plt.xlabel("Direction of Motion (degrees)")
    plt.ylabel("Firing Rate (spikes/s)")
    plt.title(title)
    
    plt.subplot(2,2,4,  polar = True)    
    spikecounts = direction_rates[:,1]
    spikecounts2 = np.append(spikecounts, direction_rates[0,1]) 
    r = np.arange(0, 361, 45)*np.pi/180
    plt.polar(r, spikecounts2,'o')
    plt.polar(fit_curve[:,0]*np.pi/180,fit_curve[:,1],'-', label="Firing Rate (spikes/s)")
    plt.title(title)
    plt.legend(loc=8)

def plot_transition_ratio(df1, df2):
    """
    plot stage transitions
    df1: normal sleep (df1 = analyse(base))
    df2: sleep depravation (df2 = analyse(depr))
    """    
    N = 5
    ind = np.arange(N)  # the x locations for the groups
    width = 0.2       # he width of the bars
    plt.close()
    plt.rc('font', family='Arial')

    fig, ax = plt.subplots(nrows=6, ncols=6, sharex='col', sharey='row')
    fig.suptitle("Comparison of the number of stage transitions (% of total transitions) (origin stage " + u'\u2192' + " dest. stage)", fontsize=20)        
    plt.subplots_adjust(wspace = 0.2,hspace = 0.4 )
    for i in range(0,6): # do not care about stage transitions > 5
        for j in range(0,6):     
            clef = '%t' + str(i) + '-' + str(j)
            normal = df1[clef].tolist()
            mean = sum(normal) / len(normal)
            normal.extend([mean])
            rects1 = ax[i,j].bar(ind, normal, width, color='b')
            depravation = df2[clef].tolist()
            mean = sum(depravation) / len(depravation)
            depravation.extend([mean])
            rects2 = ax[i,j].bar(ind+width, depravation, width, color='r')
            for label in (ax[i,j].get_xticklabels() + ax[i,j].get_yticklabels()):
                label.set_fontname('Arial')
                label.set_fontsize(8)            
            ax[i,j].set_title(str(i) + ' ' + u'\u2192' + ' ' + str(j))
            ax[i,j].set_xticks(ind+width)
            ax[i,j].set_xticklabels( ('1', '2', '3', '4', 'Avg') )
            ax[i,j].set_yticks(np.arange(0, 6, 2))
            ax[i,j].set_ylim([0,6])
    fig.legend( (rects1[0], rects2[0]), ('Baseline', 'Recovery'), loc = 'lower right', fontsize=10)

def plot_ratio_cormats():
	nprops = 19
	cormat_yn = N.loadtxt("cormat_to_do_pca.dat")
	cormat_jc = N.loadtxt("corrcoeff_jc.dat", usecols=[2])
	cormat_jc = cormat_jc.reshape(nprops, nprops)
	rat = cormat_yn/cormat_jc

	fig = plt.figure(1)
	plt.clf()
	plt.subplots_adjust(hspace=0.3, wspace=0.3)
	fontsize=8
	for ip in xrange(nprops): 
		ax = fig.add_subplot(4, 5, ip+1)
		ax.plot(N.arange(nprops), rat[ip, :], 'ro', ms=2.5)
		ax.xaxis.set_major_locator(plt_ticker.MaxNLocator(4))
		ax.yaxis.set_major_locator(plt_ticker.MaxNLocator(4))
		for tick in ax.xaxis.get_major_ticks():
			tick.label1.set_fontsize(fontsize)
		for tick in ax.yaxis.get_major_ticks():
			tick.label1.set_fontsize(fontsize)

		ii = N.where(N.abs(rat[ip, :]-1.) > 0.01)[0]
		print "-"*10
		print "{:d}: {:}".format(ip, ii)
		print cormat_yn[ip, ii]
		print cormat_jc[ip, ii]

def SVD_plot(SVStreams, SValues, stachans, title=False):
    r"""Function to plot the singular vectors from the clustering routines, one\
    plot for each stachan

    :type SVStreams: list of :class:Obspy.Stream
    :param SVStreams: See clustering.SVD_2_Stream - will assume these are\
            ordered by power, e.g. first singular vector in the first stream
    :type SValues: list of float
    :param SValues: List of the singular values corresponding to the SVStreams
    :type stachans: list
    :param stachans: List of station.channel
    """
    for stachan in stachans:
        print(stachan)
        plot_traces = [SVStream.select(station=stachan.split('.')[0],
                                       channel=stachan.split('.')[1])[0]
                       for SVStream in SVStreams]
        fig, axes = plt.subplots(len(plot_traces), 1, sharex=True)
        axes = axes.ravel()
        for i, tr in enumerate(plot_traces):
            y = tr.data
            x = np.linspace(0, len(y) * tr.stats.delta, len(y))
            axes[i].plot(x, y, 'k', linewidth=1.1)
            ylab = 'SV '+str(i+1)+'='+str(round(SValues[i] / len(SValues), 2))
            axes[i].set_ylabel(ylab, rotation=0)
            axes[i].yaxis.set_ticks([])
            print(i)
        axes[-1].set_xlabel('Time (s)')
        plt.subplots_adjust(hspace=0)
        if title:
            axes[0].set_title(title)
        else:
            axes[0].set_title(stachan)
        plt.show()
    return

    def draw(self, description=None, ofile="test.png"):

        plt.clf()

        for f in self.funcs:
            f()

        plt.axis("off")

        ax = plt.gca()
        ax.set_aspect("equal", "datalim")

        f = plt.gcf()
        f.set_size_inches(12.8, 7.2)

        if description is not None:
            plt.text(0.025, 0.05, description, transform=f.transFigure)

        if self.xlim is not None:
            plt.xlim(*self.xlim)

        if self.ylim is not None:
            plt.ylim(*self.ylim)

        plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)

        # dpi = 100 for 720p, 150 for 1080p
        plt.savefig(ofile, dpi=150)

def display_collision3D(collision):

    jets,muons,electrons,photons,met = collision

    lines = draw_beams()

    pmom = np.array(jets).transpose()[1:4].transpose()
    origin = np.zeros((len(jets),3))
    lines += draw_jet3D(origin=origin,pmom=pmom)

    pmom = np.array(muons).transpose()[1:4].transpose()
    origin = np.zeros((len(muons),3))
    lines += draw_muon3D(origin=origin,pmom=pmom)

    pmom = np.array(electrons).transpose()[1:4].transpose()
    origin = np.zeros((len(electrons),3))
    lines += draw_electron3D(origin=origin,pmom=pmom)

    pmom = np.array(photons).transpose()[1:4].transpose()
    origin = np.zeros((len(photons),3))
    lines += draw_photon3D(origin=origin,pmom=pmom)

    fig = plt.figure(figsize=(6,4),dpi=100)
    ax = fig.add_subplot(1,1,1)
    ax = fig.gca(projection='3d')
    plt.subplots_adjust(top=0.98,bottom=0.02,right=0.98,left=0.02)

    for l in lines:
        ax.add_line(l)

    ax.set_xlim(-200,200)
    ax.set_ylim(-200,200)
    ax.set_zlim(-200,200)

def trace_plot(data, ylim):
    
    clf()
    plt.subplots_adjust(left=0.15);
    
    plt.subplot(411)
    plt.plot(np.arange(0,30000*180)/30000., data[0][:30000*3*60])
    plt.xticks(alpha = 0)
    plt.yticks(fontsize = 'large')
    plt.ylabel('Voltage (mV)', size = 'x-large')
    plt.ylim(ylim)
    
    plt.subplot(412)
    plt.plot(np.arange(0,30000*180)/30000.,data[1][:30000*3*60])
    plt.yticks(alpha = 0)
    plt.xticks(alpha = 0)
    plt.ylim(ylim)
    
    plt.subplot(413)
    plt.plot(np.arange(0,30000*180)/30000.,data[2][:30000*3*60])
    plt.xticks(alpha = 0)
    plt.yticks(alpha = 0)
    plt.ylim(ylim)
    
    plt.subplot(414)
    plt.plot(np.arange(0,30000*180)/30000., data[3][:30000*3*60])
    plt.xlabel('Time (s)', size = 'x-large')
    plt.xticks(fontsize = 'large')
    plt.yticks(alpha = 0)
    plt.ylim(ylim)

	def plot_scatter(self, iclus):
		observed = rfn.read_data()
		fig = plt.figure(2)
		plt.clf()

		trueM = observed[iclus*self.nlos, 1]
		mass = observed[iclus*self.nlos:(iclus+1)*self.nlos, 2:]
		mass, masstype = self.get_valid_array_altogether(mass, trueM)

		ii = 0
		for iobs in xrange(self.nobs):
			for jobs in xrange(iobs+1, self.nobs, 1):

				ax = fig.add_subplot(4, 3, ii+1)
				ax.plot(mass[:, iobs], mass[:, jobs], 'bo', ms=2.5)

				plt.subplots_adjust(**self.adjustparam)
				ax.xaxis.set_major_locator(plt_ticker.MaxNLocator(3))
				ax.yaxis.set_major_locator(plt_ticker.MaxNLocator(3))

				#ax.set_xlabel("%s/<%s>" % (self.obsname[iobs], self.obsname[iobs])) 
				#ax.set_ylabel("%s/<%s>" % (self.obsname[jobs], self.obsname[jobs])) 

				fontsize=9
				for tick in ax.xaxis.get_major_ticks():
					tick.label1.set_fontsize(fontsize)
				for tick in ax.yaxis.get_major_ticks():
					tick.label1.set_fontsize(fontsize)

				ii += 1
				print self.obsname[iobs], self.obsname[jobs]

		plt.savefig(os.path.join("paper", "figure", "scatter_%s_clus%d.eps" % (masstype, iclus)), orientation='portrait', transparent=True)

    def run(self):
        lines = open(self.inFilename).readlines()
        data = []

        for line in lines:
            data.append(float(line.strip()))
        x = np.asarray(data)
        fig = plt.figure(figsize=(7, 3))        
     
        ax = fig.add_subplot(111)
        plt.subplots_adjust(left = 0.15, bottom = 0.15, wspace = 0)
        plt.xlabel(self.options.xlab)
       
        
        plt.ylabel(self.options.ylab)
        if self.options.logy == True:
            ax.set_yscale('log')
        plt.title(self.options.title)
        if self.options.plotType == 'hist':
            plt.xlim(0,x.max())
            n,bins,patches = plt.hist(x, self.options.bins,  histtype='bar',
                                 color=['crimson'],normed=False, alpha=0.85)
        else:
            plt.xlim(0,x.size)
            line, = plt.plot(range(x.size), x, 'r-', label = self.options.label)
        if self.options.label:
            ax.legend()
        plt.savefig(self.outFilename)

def plot_spectrograms(data, rate, subject, condition):
    """
    Creates spectrogram subplots for all 9 channels
    """
    fig = plt.figure()   
   # common title
    fname = 'Spectrogram - '+'Subject #'+subject+' '+condition+' Dataset'
    fig.suptitle(fname, fontsize=14, fontweight='bold')            
    # common ylabel
    fig.text(0.06, 0.5, 'ylabel', 
             ha='center', va='center', rotation='vertical',
             fontsize=14, fontweight='bold')
    # use this to stack EEG, EOG, EMG on top of each other         
    sub_order = [1,4,7,10,2,5,3,6,9]          
    
    for ch in range(0, len(data)):
        plt.subplot(4, 3, sub_order[ch])
        plt.subplots_adjust(hspace=.6)  # adds space between subplots
        plt.title(channel_name[ch])    
        Pxx, freqs, bins, im = plt.specgram(data[ch],NFFT=512,Fs=rate)       
        plt.ylim(0,70)
        plt.xlabel('Time (Seconds)')
        plt.ylabel('Frequency (Hz)')
    
    #fig.savefig(fname+'.pdf', format='pdf')  buggy resolution problem
    return

def plot_proj_scatter(pccoefs, eval, figname):
	"""
	plot scatters using the values projected on pc
	"""
	fig = plt.figure(4)
	plt.clf()

	eval /= N.sum(eval)

	param = dict(fontsize='small')
	
	range = N.max(N.abs(pccoefs))
	labelx = -0.15

	nplt = 4
	ncol = 3
	nrow = 2
	iplt = 1
	for iprops in xrange(nplt):
		for jprops in xrange(iprops+1, nplt):
			ax = fig.add_subplot(nrow, ncol, iplt)
			plt.subplots_adjust(hspace=0.2, wspace=0.25, left=0.08, right=0.95, bottom=0.08, top=0.9)
			ax.plot(pccoefs[:, iprops], pccoefs[:, jprops], 'ro', ms=3.)
			plt.axis('equal')
			#ax.set_xlim(-range, range)
			#ax.set_ylim(-range, range)
			ax.set_xlabel('PC%d(%.2f)' % (iprops, eval[iprops]), **param)
			ax.set_ylabel('PC%d(%.2f)' % (jprops, eval[jprops]), **param)
			ax.yaxis.set_label_coords(labelx, 0.5)
			iplt += 1

	plt.savefig(figname)

def plot_transition(df1, df2):
    """
    plot stage transitions
    df1: normal sleep (df1 = analyse(base))
    df2: sleep depravation (df2 = analyse(depr))
    """    
    N = 5
    ind = np.arange(N)  # the x locations for the groups
    width = 0.2       # the width of the bars
    plt.close()

    fig, ax = plt.subplots(nrows=6, ncols=6, sharex='col', sharey='row')
    plt.subplots_adjust(wspace = 0.2,hspace = 0.4 )
    for i in range(0,6): # do not care about stage transitions > 5
        for j in range(0,6):     
            clef = 't' + str(i) + '-' + str(j)
            normal = df1[clef].tolist()
            mean = sum(normal) / len(normal)
            normal.extend([mean])
            rects1 = ax[i,j].bar(ind, normal, width, color='r')
            depravation = df2[clef].tolist()
            mean = sum(depravation) / len(depravation)
            depravation.extend([mean])
            rects2 = ax[i,j].bar(ind+width, depravation, width, color='y')
            ax[i,j].set_title('t' + str(i) + '-' + str(j))
            ax[i,j].set_xticks(ind+width)
            ax[i,j].set_xticklabels( ('1', '2', '3', '4', 'Avg') )

def paired_boxplot_o(boxes):
    """
    Wrapper around plt.boxplot to draw paired boxplots
    for a set of boxes. 
    
    Input is the same as plt.boxplot:
        Array or a sequence of vectors.
    """
    fig = plt.figure(figsize=(len(boxes) / 2.5, 4))
    ax1 = fig.add_subplot(111)
    plt.subplots_adjust(left=0.075, right=0.95, top=0.9, bottom=0.25)
    bp = ax1.boxplot(boxes, notch=0, positions=np.arange(len(boxes)) + 
                     1.5 * (np.arange(len(boxes)) / 2), patch_artist=True)
    [p.set_color(colors[0]) for p in bp['boxes'][::2]]
    [p.set_color('black') for p in bp['whiskers']]
    [p.set_color('black') for p in bp['fliers']]
    [p.set_alpha(.4) for p in bp['fliers']]
    [p.set_alpha(.6) for p in bp['boxes']]
    [p.set_edgecolor('black') for p in bp['boxes']]
    ax1.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
                  alpha=0.5)
    
    # Hide these grid behind plot objects
    ax1.set_axisbelow(True)
    ax1.set_ylabel('$Log_{2}$ RNA Expression')
    ax1.set_xticks(3.5 * np.arange(len(boxes) / 2) + .5)
    return ax1, bp

 def segmentation(self, threshold):
     img = self.spectrogram["data"]
     mask = (img > threshold).astype(np.float)
     hist, bin_edges = np.histogram(img, bins=60)
     bin_centers = 0.5*(bin_edges[:-1] + bin_edges[1:])
     binary_img = mask > 0.5
     plt.figure(figsize=(11,8))
     plt.subplot(131)
     plt.imshow(img)
     plt.axis('off')
     plt.subplot(132)
     plt.plot(bin_centers, hist, lw=2)
     print(threshold)
     plt.axvline(threshold, color='r', ls='--', lw=2)
     plt.text(0.57, 0.8, 'histogram', fontsize=20, transform = plt.gca().transAxes)
     plt.text(0.45, 0.75, 'threshold = '+ str(threshold)[0:5], fontsize=15, transform = plt.gca().transAxes)
     plt.yticks([])
     plt.subplot(133)     
     plt.imshow(binary_img)
     plt.axis('off')
     plt.subplots_adjust(wspace=0.02, hspace=0.3, top=1, bottom=0.1, left=0, right=1)
     plt.show()
     print(img.max())
     print(binary_img.max())
     
     return mask

def animate_plot(w_dir=None, **kwargs):
    if w_dir == None: w_dir=os.getcwd()

    cdir = os.getcwd()
    files = []
    os.chdir(w_dir)
    os.system('mkdir j_movie')
    D0 = plp.pload(0,w_dir=w_dir)
    fig = plt.figure(num=1,figsize=[7,7])
    ax = fig.add_subplot(111)
    plt.subplots_adjust(left=0.1, bottom=0.15)
    frnum = kwargs.get('frames',plp.time_info(w_dir=w_dir)['time'])
    
    
    
    for i in [0,10,25,50,150,250,319,325,340,450,638]:  # 50 frames
        D = plp.pload(i,w_dir=w_dir)
        Ra = da.Rad_Average()
        xitem = D.x1
	yitem = D.x1*D.v3[:,32]
    #    yitem = Ra.Sigma(D,ul=1.0,urho=1.0e-9,Mstar=10.0,Gammae=5.0/3.0)
        ax.cla()
        
        
        

        if kwargs.get('pltype','normal') == 'normal':
            ax.plot(xitem,yitem,'k-')
        if kwargs.get('pltype','normal') == 'logx':
            ax.semilogx(xitem,yitem,'k-')
        if kwargs.get('pltype','normal') == 'logy':
            ax.semilogy(xitem,yitem,'k-')
        if kwargs.get('pltype','normal') == 'loglog':
            ax.loglog(xitem,yitem,'k-')

        ax.minorticks_on()
        ax.set_xlabel(kwargs.get('xlabel',r'XLabel'))
        ax.set_ylabel(kwargs.get('ylabel',r'YLabel'))
        ax.set_title(kwargs.get('title',r'Title'))
        
        
        #ax.text(90.0,1.3,r'$N_{\rm rot} = %04d$'%i)
        axcolor = 'white'
        axtime = plt.axes([0.1, 0.01, 0.8, 0.02], axisbg=axcolor)
        axtime.cla()
        stime = Slider(axtime, 'Nsteps', 0, frnum , valinit=D.time)

        #plt.draw()
        fname = w_dir+'j_movie/_tmp%03d.png'%i
        
        print 'Saving frame', fname
        fig.savefig(fname)
        files.append(fname)
        
    

    print 'Making movie animation.mpg - this make take a while'
    os.system("mencoder 'mf://j_movie/_tmp*.png' -mf type=png:fps=10 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o j_movie/animation.mpg")
    os.chdir(cdir)