def plotDataAndFit_S():
    """
	
		Plot \\xi( |s| )
	
	"""

    binSize = 4.0

    sss = numpy.arange(0.0, 50.0) * 4.0 + 2.0
    yyy_dat = numpy.zeros(shape=(50, 2))
    yyy_fit = numpy.zeros(shape=(50, 2))

    cut = numpy.logical_and((xxx__ != 0.0), (xxx__ < 200.0))
    xxx = xxx__[cut]
    xi_dat = xi_dat__[cut]
    xi_fit = xi_fit__[cut]
    xi_err = 1.0 / (xi_err__[cut] ** 2.0)

    for i in range(0, xxx.size):
        sIdx = int(xxx[i] / binSize)
        yyy_dat[sIdx][0] += xi_dat[i] * xi_err[i]
        yyy_dat[sIdx][1] += xi_err[i]
        yyy_fit[sIdx][0] += xi_fit[i] * xi_err[i]
        yyy_fit[sIdx][1] += xi_err[i]

    yyy_dat[:, 0] /= yyy_dat[:, 1]
    yyy_dat[:, 1] = numpy.sqrt(1.0 / yyy_dat[:, 1])
    yyy_fit[:, 0] /= yyy_fit[:, 1]
    yyy_fit[:, 1] = numpy.sqrt(1.0 / yyy_fit[:, 1])

    ### Plot the results

    for i in numpy.arange(0, 3):

        a = ""
        if i == 1:
            a += "|s|."
        elif i == 2:
            a += "|s|^{2}."

        coef = numpy.power(sss, 1.0 * i)

        plt.errorbar(
            sss,
            coef * yyy_dat[:, 0],
            yerr=coef * yyy_dat[:, 1],
            linestyle="",
            marker="o",
            color="blue",
            label=r"$<Simu>$",
        )
        plt.errorbar(sss, coef * yyy_fit[:, 0], color="red", label=r"$<Fit>$")
        plt.xlabel(r"$|s| \, [h^{-1} Mpc]$")
        plt.ylabel(r"$" + a + "\\xi(|s|)$")
        myTools.deal_with_plot(False, False, True)

        plt.show()

    return

def PlotPatchBins(Sc, PatchData, NumBins, color, MinimumBinSize=7, ErrorBars=True):
    """
    Plot E*R* data binned from the hilltop pacth data.
    """
    E_s = E_Star(Sc, PatchData[6], PatchData[2])
    R_s = R_Star(Sc, PatchData[10], PatchData[2])

    bin_x, bin_std_x, bin_y, bin_std_y, std_err_x, std_err_y, count = Bin.bin_data_log10(E_s, R_s, NumBins)

    # filter bins based on the number of data points used in their calculation
    bin_x = np.ma.masked_where(count < MinimumBinSize, bin_x)
    bin_y = np.ma.masked_where(count < MinimumBinSize, bin_y)
    # these lines produce a meaningless warning - don't know how to solve it yet.

    if ErrorBars:
        # only plot errorbars for y as std dev of x is just the bin width == meaningless
        plt.scatter(
            bin_x,
            bin_y,
            c=count,
            s=50,
            edgecolor="",
            cmap=plt.get_cmap("autumn_r"),
            label="Binned Patch Data",
            zorder=100,
        )
        plt.errorbar(bin_x, bin_y, yerr=std_err_y, fmt=None, ecolor="k", elinewidth=2, capsize=3, zorder=0)
        cbar = plt.colorbar()
        cbar.set_label("Number of values per bin")

    else:
        plt.errorbar(bin_x, bin_y, fmt="o", color=color, label="No. of Bins = " + str(NumBins))

def display(params_estimated):

    # Construct matrix of experimental data and variance columns of interest
    exp_obs_norm = exp_data[data_names].view(float).reshape(len(exp_data), -1).T
    var_norm = exp_data[var_names].view(float).reshape(len(exp_data), -1).T
    std_norm = var_norm ** 0.5

    # Simulate model with new parameters and construct a matrix of the
    # trajectories of the observables of interest, normalized to 0-1.
    solver.run(params_estimated)
    obs_names_disp = obs_names + ['aSmac']
    sim_obs = solver.yobs[obs_names_disp].view(float).reshape(len(solver.yobs), -1)
    totals = obs_totals + [momp_obs_total]
    sim_obs_norm = (sim_obs / totals).T

    # Plot experimental data and simulation on the same axes
    colors = ('r', 'b')
    for exp, exp_err, sim, c in zip(exp_obs_norm, std_norm, sim_obs_norm, colors):
        plt.plot(exp_data['Time'], exp, color=c, marker='.', linestyle=':')
        plt.errorbar(exp_data['Time'], exp, yerr=exp_err, ecolor=c,
                     elinewidth=0.5, capsize=0, fmt=None)
        plt.plot(solver.tspan, sim, color=c)
    plt.plot(solver.tspan, sim_obs_norm[2], color='g')
    plt.vlines(momp_data[0], -0.05, 1.05, color='g', linestyle=':')
    plt.show()

def build_plot(profilerResults):
    # Calculate each value.
    x = []
    mean = []
    std = []
    for t in xrange(profilerResults.getLookBack()*-1, profilerResults.getLookForward()+1):
        x.append(t)
        values = np.asarray(profilerResults.getValues(t))
        mean.append(values.mean())
        std.append(values.std())

    # Cleanup
    plt.clf()
    # Plot a line with the mean cumulative returns.
    plt.plot(x, mean, color='#0000FF')

    # Error bars starting on the first lookforward period.
    lookBack = profilerResults.getLookBack()
    firstLookForward = lookBack+1
    plt.errorbar(
        x=x[firstLookForward:], y=mean[firstLookForward:], yerr=std[firstLookForward:],
        capsize=3,
        ecolor='#AAAAFF', alpha=0.5
    )

    # Horizontal line at the level of the first cumulative return.
    plt.axhline(
        y=mean[lookBack],
        xmin=-1*profilerResults.getLookBack(), xmax=profilerResults.getLookForward(),
        color='#000000'
    )

    plt.xlim(profilerResults.getLookBack()*-1-0.5, profilerResults.getLookForward()+0.5)
    plt.xlabel('Time')
    plt.ylabel('Cumulative returns')

def plot_categorical_scatter_with_mean(vals, categoryLabels, jitter=True, colours=None, xlabel=None, ylabel=None, title=None):
    import matplotlib.colors
    import scipy.stats
    import pdb
    numCategories = len(vals)
    plt.hold(True)
    if colours is None:
        colours = plt.cm.gist_rainbow(np.linspace(0,1,numCategories))
    for category in range(numCategories):
        edgeColour = matplotlib.colors.colorConverter.to_rgba(colours[category], alpha=0.5)
        xval = (category+1)*np.ones(len(vals[category]))
        if jitter:
            jitterAmt = np.random.random(len(xval))
            xval = xval + (0.3 * jitterAmt) - 0.15
        #pdb.set_trace()
        plt.plot(xval, vals[category], 'o', mec=edgeColour, mew = 4, mfc='none', ms=16)
        mean = np.mean(vals[category])
        sem = scipy.stats.sem(vals[category])
        print mean, sem
        plt.plot(category+1, mean, 'o', color='k', mec=colours[category], ms=20)
        plt.errorbar(category+1, mean, yerr = sem, color=colours[category])
    plt.xlim(0,numCategories+1)
    plt.ylim(0,1)
    ax = plt.gca()
    ax.set_xticks(range(1,numCategories+1))
    ax.set_xticklabels(categoryLabels, fontsize=16)
    if xlabel is not None:
        plt.xlabel(xlabel, fontsize=20)
    if ylabel is not None:
        plt.ylabel(ylabel, fontsize=20)
    if title is not None:
        plt.title(title)
    plt.show()

	def fitplot(self,pars):
		
		AnalyticTTVs = []
		for parset in pars:
			m,m1,ex,ey,ex1,ey1 = parset
			AnalyticTTVs.append( self.get_ttvs(ex,ey,ex1,ey1) )
#		#
		pl0tr = self.transits
		pl1tr = self.transits1
		N = self.trN
		N1 = self.trN1
		errs,errs1 = self.input_data[:,2],self.input_data1[:,2]
#		#
		symbols = ['x','o','d']
		## Figure 1 ##
		plt.figure()
		plt.subplot(211)
		plt.errorbar(pl0tr, pl0tr - self.p*N - self.T0,yerr=errs,fmt='ks')
		for i,ttvs in enumerate(AnalyticTTVs):
			plt.plot(pl0tr , ttvs[0] * m1 ,'k%s'% symbols[i%len(symbols)] ) 
		plt.subplot(212)
		plt.errorbar(pl1tr , pl1tr - self.p1*N1 - self.T10 ,yerr=errs1,fmt='rs')
		for i,ttvs in enumerate(AnalyticTTVs):
			plt.plot(pl1tr , ttvs[1] * m  ,'r%s'% symbols[i%len(symbols)] ) 
		plt.show()

def plotres(psr,deleted=False,group=None,**kwargs):
    """Plot residuals, compute unweighted rms residual."""

    res, t, errs = psr.residuals(), psr.toas(), psr.toaerrs
    
    if (not deleted) and N.any(psr.deleted != 0):
        res, t, errs = res[psr.deleted == 0], t[psr.deleted == 0], errs[psr.deleted == 0]
        print("Plotting {0}/{1} nondeleted points.".format(len(res),psr.nobs))

    meanres = math.sqrt(N.mean(res**2)) / 1e-6
    
    if group is None:
        i = N.argsort(t)
        P.errorbar(t[i],res[i]/1e-6,yerr=errs[i],fmt='x',**kwargs)
    else:
        if (not deleted) and N.any(psr.deleted):
            flagmask = psr.flagvals(group)[~psr.deleted]
        else:
            flagmask = psr.flagvals(group)

        unique = list(set(flagmask))
            
        for flagval in unique:
            f = (flagmask == flagval)
            flagres, flagt, flagerrs = res[f], t[f], errs[f]
            i = N.argsort(flagt)
            P.errorbar(flagt[i],flagres[i]/1e-6,yerr=flagerrs[i],fmt='x',**kwargs)
        
        P.legend(unique,numpoints=1,bbox_to_anchor=(1.1,1.1))

    P.xlabel('MJD'); P.ylabel('res [us]')
    P.title("{0} - rms res = {1:.2f} us".format(psr.name,meanres))

def draw_zprofile(
        analyzed_root_files,
        energy,
        histlabel=None,
        histcolor = '#000000'):
    
    nmodules = 14
    nsectors = 16
    cut_event_numb = 0. # counter for number of events after cuts are applied

    # add to tree all edm_analyzed reco files for electrongun
    rh_tree = ROOT.TChain("demo/rh_tree")
    rh_tree.Add(analyzed_root_files)

    # initialize arrays for module energies
    total_module_energy = nmodules * [0.] # array with energy of all sectors for each module, air
    total_module_energy_per_event = nmodules * [0.]

    # initialize lists for error calculation
    # sigma = sqrt(<x^2> - <x>^2), mean_error = sigma/sqrt(2)
    module_energies_squared = nmodules * [0.] # sum of squares of module energies 
    bin_error =  nmodules * [0.]


    # loop over all events 

    for i in range(rh_tree.GetEntries()):
        rh_tree.GetEntry(i)
        px, py, pz = rh_tree.gen_part_momentum_x, rh_tree.gen_part_momentum_y, rh_tree.gen_part_momentum_z
        momentum_vector = ROOT.TVector3(px, py, pz)
        eta = momentum_vector.Eta()
        if (eta > -6.4) & (eta < -5.5):
            cut_event_numb += 1
            for module_numb in range(nmodules):
                total_module_energy_per_event[module_numb] = 0
                for sector_numb in range(nsectors):
                    total_module_energy[module_numb] += rh_tree.energy_castor[sector_numb*nmodules + module_numb]
                    total_module_energy_per_event[module_numb] += rh_tree.energy_castor[sector_numb*nmodules + module_numb]
                module_energies_squared[module_numb] += total_module_energy_per_event[module_numb]**2
            
    print "cut event number: ", cut_event_numb


    for module_numb in range(nmodules):
        bin_error[module_numb] = 1./energy*1./np.sqrt(cut_event_numb) * np.sqrt(module_energies_squared[module_numb]/cut_event_numb
                                - (total_module_energy[module_numb]/cut_event_numb)**2)
        print "module; ", module_numb+1, ", bin_error: " , bin_error[module_numb]

    mean_module_energy = np.array(total_module_energy)/cut_event_numb


    # draw histogram 
    # plt.bar(np.arange(1,15,1), mean_module_energy/energy,
    #         yerr=np.array(bin_error)/energy, # value of error bars
    #         width=1, fill=False, edgecolor=histedgecolor, alpha = transparency, label=str(energy)+" GeV",
    #         error_kw=dict(ecolor=histedgecolor)) # set collor of error bar
    plt.hist(np.arange(1.,15.,1), bins=np.arange(1,16), weights=mean_module_energy/energy,
             histtype=u'step', align = u'mid', color=histcolor , label=histlabel)
    plt.errorbar(np.arange(1.5,15.5,1), mean_module_energy/energy, yerr=np.array(bin_error),
                 fmt='none', ecolor=histcolor)

	def plot_minos_fit(self,p,decay="X",title="Fit Results",erange=9,step=4.07,lum=4200):
		fig = plt.figure(figsize=(8,6))
		plt.errorbar(self.x,self.y,self.yerr,fmt='o')
		M = p[0][0]
		G = p[1][0]
		B = p[2][0]
		dMl = p[0][1]
		dMu = p[0][2]
		dGl = p[1][1]
		dGu = p[1][2]
		dBl = p[2][1]
		dBu = p[2][2]
		x_fit = np.linspace(min(self.x),max(self.x),num=100)
		plt.plot(x_fit,self.convBWG(x_fit,M,G,B))
		plt.xlabel("$\sqrt{\hat{s}} (MeV)$",fontsize=16)
		plt.ticklabel_format(useOffset=False)
		plt.ylabel("Counts",fontsize=16)
		plt.title(title,fontsize=16)
		lbl1 = "Input:\n$\mathcal{L}=%d pb^{-1}$\n$\Delta=%.3f\ MeV$\n$\delta\sqrt{\hat{s}} = %.3f MeV$" % (lum,step,self.beam)
		lbl1 = lbl1 + "\n$M_h = 125.0 GeV$\n$\Gamma_h = 4.07 MeV$\n$Br(h^0\\rightarrow$%s$) = %.3f$" % (decay, self.higgs[2])
		lbl1 = lbl1 + "\n$\sigma_{bkg} = %.2f pb^{-1}$" % (self.bkg)
		lbl2 = "\nFit results:\n"
		lbl2 = lbl2 + "$\Delta M_h = %.3f_{-%.3f}^{+%.3f}\ MeV$\n" % (M-self.higgs[0], -1*dMl, dMu)
		lbl2 = lbl2 + "$\Gamma_h = %.3f_{-%.3f}^{+%.3f} \ MeV$\n" % (G, -1*dGl, dGu)
		lbl2 = lbl2 + "$Br(h^0\\rightarrow$%s$) = %.3f_{-%.3f}^{+%.3f}$\n" % (decay, B, -1*dBl, dBu)
		plt.annotate(lbl1, [0.1,0.6], xycoords='axes fraction',fontsize=15)
		plt.annotate(lbl2, [0.7,0.6], xycoords='axes fraction',fontsize=15)
		return plt

def plot_ts(table):
    """
    Plot dh and dAGC time series and the correlation dAGC x dh.
    """
    sys.path.append('/Users/fpaolo/code/misc')
    from util import poly_fit
    # load data from Table
    time2 = table.cols.time2[:] 
    month = table.cols.month[:] 
    dh_mean = table.cols.dh_mean[:] 
    dh_error = table.cols.dh_error[:] 
    dg_mean = table.cols.dg_mean[:] 
    dg_error = table.cols.dg_error[:] 
    dates = [dt.datetime(y, m, 15) for y, m in zip(time2, month)]
    # plot TS
    fig = plt.figure()
    plt.subplot(211)
    plt.errorbar(dates, dh_mean, yerr=dh_error, linewidth=2)
    plt.ylabel('dh (m)')
    plt.subplot(212)
    plt.errorbar(dates, dg_mean, yerr=dg_error, linewidth=2)
    plt.ylabel('dAGC (dB)')
    fig.autofmt_xdate()
    # plot correlation
    dg_fit, dh_fit, _ = poly_fit(dg_mean, dh_mean)
    plt.figure()
    plt.plot(dg_mean, dh_mean, 'o')
    plt.plot(dg_fit, dh_fit, linewidth=2.5)
    plt.xlabel('dAGC (dB)')
    plt.ylabel('dh (m)')
    corr = np.corrcoef(dg_mean, dh_mean)[0,1]
    print 'correlation = %.2f' % corr

def mw_cdf(x_hist_vals, hist_vals, a_coeff, figs, plot=False):
    max_a = np.sum(hist_vals)
    area_a = np.ones(len(hist_vals))
    for el in xrange(len(hist_vals)):
        area_a[el] = np.sum(hist_vals[0:el+1])

    
    c_d_f = area_a/max_a
    interp = interp1d(c_d_f, x_hist_vals)
    a_best = interp(0.5)
    a_limits = interp(np.array([0.5 - 0.683/2.0, 0.5 + 0.683/2.0]))

    decim = [math.trunc(np.abs(np.log10(a_best - a_limits[0])))+2, math.trunc(np.abs(np.log10(a_limits[1] - a_best)))+2]
             
    uncertainties = np.array([round(a_best - a_limits[0], decim[0]), round(a_limits[1] - a_best, decim[1])])

    if plot:
        plt.figure(figs)
        figs += 1
        plt.clf()
        plt.scatter(x_hist_vals, c_d_f, marker='+')
        plt.plot((a_best, a_best), (( c_d_f.max(), 0)), 'g')
        plt.errorbar(a_best, 0.5, xerr=[[uncertainties[0]], [uncertainties[1]]], fmt='^', color='red')
        plt.ylabel('CDF ')
        plt.xlabel(r'a$_'+str(a_coeff)+'$ values')
        
        plt.title(r'Result: a$_'+str(a_coeff)+' = '+str(round(a_best, np.max(decim)))+'_{-'+str(uncertainties[1])+'}^{+'+str(uncertainties[0])+'}$')
        plt.show() #in most cases unnecessary
        
    return figs

def plot_cluster_size(file, outfile='cluster_size.pdf', title=None, grid=False, epochs=[],
                      labels=[], errorbars=True):
    reader = csv.reader(file)
    rownum = 0
    data_read = False

    for row in reader:
        # Skip commented lines
        if re.match('^\s*#', row[0]) != None:
            continue

        rownum += 1

        if rownum == 1:
            colnames = row[1:]

            for i in range(len(colnames)):
                colnames[i] = re.sub('_', ' ', colnames[i])

            continue
        else:
            row = list(map(float, row))

            if (len(epochs) == 0 or row[0] in epochs):
                if data_read:
                    data = numpy.vstack((data,row))
                else:
                    data = numpy.array(row)
                    data_read = True

    if data_read:
        fig = plt.figure()

        if len(labels) > 0 and len(labels) == len(colnames):
            colnames = labels

        plot_cols = list(range(2, data.shape[1], 3))

        # Plot the number of clusters
        for t in plot_cols:
            if errorbars:
                e = data[:, t+1]
            else:
                e = None

            plt.errorbar(data[:,0], data[:,t], yerr=e, xerr=None, label=string.capitalize(colnames[t-1]))

        if grid:
            plt.grid()

        plt.xlabel("Time (epoch)")
        plt.ylabel("Cluster Size (cells)")

        if title:
            plt.title(title)

        plt.legend(loc=0)
        plt.savefig(outfile)
    else:
        print("Could not generate plot: No data match given parameters")

def histograma(h,l):

    Hist = genhistograma(h)
    
    if h <= 3:
        x = np.arange(1.5,13,1)
        H = np.histogram(Hist, bins = [1,2,3,4,5,6,7,8,9,10,11,12,13])
        y = H[0]
        plt.hist(Hist, bins = [1,2,3,4,5,6,7,8,9,10,11,12,13])
        plt.grid(True)
        plt.xlim((1,13))
        plt.ylabel("frecuencia")
        plt.title("Histograma %d"%(h+1))
        plt.errorbar(x,y, yerr = np.sqrt(y), fmt = '.')
        plt.savefig("hist%d"%(h))
        plt.close()
        
    elif h > 3:
        x = np.arange(0.5,11,1)
        H = np.histogram(Hist, bins = [0,1,2,3,4,5,6,7,8,9,10,11])
        y = H[0]
        plt.hist(Hist, bins = [0,1,2,3,4,5,6,7,8,9,10,11])
        plt.grid(True)
        plt.xlim((0,11))
        plt.ylabel("frecuencia")
        plt.title("Histograma %d"%(h+1))
        plt.errorbar(x,y, yerr = np.sqrt(y), fmt = '.')
        plt.savefig("hist%d"%(h))
        plt.close()

def plot_k2_curve(k2_arr, conc_list):
    k2_means = np.mean(k2_arr, axis=0)
    k2_sds = np.std(k2_arr, axis=0)
    # Plot k2 on a linear scale
    plt.figure()
    plt.errorbar(conc_list, k2_means, yerr=k2_sds / np.sqrt(3))
    plt.title('$k_2$')

def fit_and_plot_gains(gains, label, label_ypos, color, show_data=True, effective=True, boost=None, final=False):
    """Fit and plot a bunch of gains."""
    import matplotlib.pyplot as plt
    
    
    if final:
        expected_gain = np.asarray(gains['effective gain'] / 2.7)
    elif effective:
        expected_gain = np.asarray(gains['effective gain'])
    else:
        expected_gain = np.asarray(gains['gain'])
    model_gain = np.asarray(gains['fit gain'])
    model_noise = np.asarray(gains['fit sigma'])
    
    y = model_gain * np.sqrt(1.0/model_noise)
    A = np.vstack([expected_gain, np.zeros(len(expected_gain))]).T
    A *= np.sqrt(1.0/model_noise)[:,None]
    print(A, y)
    m, c = np.linalg.lstsq(A, y)[0]
    
    if show_data:
        plt.errorbar(expected_gain, model_gain, yerr=model_noise, fmt='.', ls='none', label=label, color=color)
    
    x = np.linspace(0.0, 2.0, 50)
    plt.plot(x, x * m + c, '-', label="{} Fit: $m={:.2f}$ $c={:.2f}$".format(label, m, c), color=color, alpha=0.3)
    
    if boost is not None:
        eboost = float(boost) / float(m)
        plt.text(0.98, 0.98, "fit boost: {:.1f}".format(eboost), transform=plt.gca().transAxes, ha='right', va='top')
    return m, c

def main():
    if not len(sys.argv) == 2:
        print "Usage: ./generate_plots.py estimator_name"
        exit(1)

    estimator_name = sys.argv[1]
    f = open("data/" + estimator_name + ".txt", 'r')

    x = []
    y = []
    e = []

    for line in f:
        values = line.split()
        x.append(float(values[0]))
        y.append(float(values[1]))
        e.append(1.96 * float(values[2]))

    f.close()

    plt.figure()
    plt.title(estimator_name)
    plt.xlabel('N')
    plt.ylabel('N_hat')
    plt.axis([1, 5000, 1, 7000])
    plt.plot(x, x, 'r-', label="Actual")
    plt.legend()
    plt.errorbar(x, y, yerr=e, fmt='ko')
    #plt.show()
    plt.savefig("figures/" + estimator_name + "_plot.png", format="png")

def PlotEDepSummary(gFiles,nFiles,figureName='EDepSummary.png',tParse=GetThickness,
  histKey='eDepHist'):
  """ PlotEDepSummary
  Plotss the energy deposition summary
  """
  # Extrating the average values
  gT = list()
  gDep = list()
  gDepError = list()
  nT = list()
  nDep = list()
  nDepError = list()
  for fname in gFiles:
    f = TFile(fname,'r')
    hist = f.Get(histKey)
    gT.append(GetThickness(fname))
    gDep.append(hist.GetMean())
    gDepError.append(hist.GetMeanError())
  for fname in nFiles:
    f = TFile(fname,'r')
    hist = f.Get(histKey)
    nT.append(GetThickness(fname))
    nDep.append(hist.GetMean())
    nDepError.append(hist.GetMeanError())
  # Plotting
  plt.errorbar(gT,gDep,yerr=gDepError,fmt='r+')
  plt.hold(True)
  plt.errorbar(nT,nDep,yerr=nDepError,fmt='go')
  plt.xlabel("Thickness (mm)")
  plt.ylabel("Average Energy Deposition (MeV)")
  plt.legend(["Co-60","Cf-252"])
  plt.xscale("log")
  plt.yscale("log")
  plt.grid(True)
  plt.savefig(figureName)

def multi_trial_stats():
    """
    Note: a lot of this should probably be in helper function tbh
    May want to add recursive copy of trial folders into summary folder, idk
    I also forgot to take sensitivity into account, I'll fix that once this all works
    """
    # Ask user for directory with all data directories to be used inside
    selected_dir = tkFileDialog.askdirectory()
    sub_dirs = [os.path.join(selected_dir, name) for name in os.listdir(selected_dir) if os.path.isdir(os.path.join(selected_dir, name))]

    # Create a TrialSummary Object for each trial
    trial_summaries = []
    for trial_dir in sub_dirs:
        trial_summaries.append(TrialSummary(trial_dir))

    # Create Plot of Averages and save
    plt.figure(2)  # off screen figure
    plt.clf()
    plt.grid(True)
    for summary in trial_summaries:
        x, y = summary.get_average_response()
        plt.plot(x, y, label=summary.name)
    plt.legend()
    plt.savefig(os.path.join(selected_dir, 'average_plot.png'))

    # Check if each trial can be given a number based on sample name
    summary_plot_numbered = []
    flag = True
    for trial in trial_summaries:
        try:
            num = float(trial.name)
            summary_plot_numbered.append((num, trial))
        except:
            flag = False
            break

    if flag:
        # Sort list by associated trial number (lest to greatest)
        summary_plot_numbered = sorted(summary_plot_numbered, key=lambda i: i[0])

    # Create the Summary Plots
    measurement_title = {0: "Amplitude (V)", 2: "Max amplitude current (A)", 4: "Min amplitude current (A)", 6: "Peak distance (V)", 8: "Peak separation (A)"}
    for measurement in range(0, 9, 2):
        plt.figure(2)
        plt.clf()
        values = [data.get_stats_list()[measurement] for number, data in summary_plot_numbered]
        errors = [data.get_stats_list()[measurement + 1] for number, data in summary_plot_numbered]
        plt.errorbar(range(len(values)), values, yerr=errors)
        for (label, x, y) in zip([x[0] for x in summary_plot_numbered], range(len(values)), values):
            plt.annotate(label, xy=(x, y), xytext=(5, 5), textcoords="offset points")
        plt.title(measurement_title[measurement])
        plt.xlabel("Sample")
        plt.ylabel("Value and error (one standard deviation)")
        plt.xlim([-1, len(values)])
        plt.savefig(os.path.join(selected_dir, measurement_title[measurement]))

    # Create Summary.txt
    for number, trial in summary_plot_numbered:
        trial.save_data(selected_dir)
    return

def test_complete():
    fig = plt.figure('Figure with a label?', figsize=(10, 6))

    plt.suptitle('Can you fit any more in a figure?')

    # make some arbitrary data
    x, y = np.arange(8), np.arange(10)
    data = u = v = np.linspace(0, 10, 80).reshape(10, 8)
    v = np.sin(v * -0.6)

    plt.subplot(3, 3, 1)
    plt.plot(list(xrange(10)))

    plt.subplot(3, 3, 2)
    plt.contourf(data, hatches=['//', 'ooo'])
    plt.colorbar()

    plt.subplot(3, 3, 3)
    plt.pcolormesh(data)

    plt.subplot(3, 3, 4)
    plt.imshow(data)

    plt.subplot(3, 3, 5)
    plt.pcolor(data)

    plt.subplot(3, 3, 6)
    plt.streamplot(x, y, u, v)

    plt.subplot(3, 3, 7)
    plt.quiver(x, y, u, v)

    plt.subplot(3, 3, 8)
    plt.scatter(x, x**2, label='$x^2$')
    plt.legend(loc='upper left')

    plt.subplot(3, 3, 9)
    plt.errorbar(x, x * -0.5, xerr=0.2, yerr=0.4)

    ###### plotting is done, now test its pickle-ability #########

    # Uncomment to debug any unpicklable objects. This is slow (~200 seconds).
#    recursive_pickle(fig)

    result_fh = BytesIO()
    pickle.dump(fig, result_fh, pickle.HIGHEST_PROTOCOL)

    plt.close('all')

    # make doubly sure that there are no figures left
    assert_equal(plt._pylab_helpers.Gcf.figs, {})

    # wind back the fh and load in the figure
    result_fh.seek(0)
    fig = pickle.load(result_fh)

    # make sure there is now a figure manager
    assert_not_equal(plt._pylab_helpers.Gcf.figs, {})

    assert_equal(fig.get_label(), 'Figure with a label?')

def plot_magnitude_time_scatter(catalogue, plot_error=False, filename=None,
        filetype='png', dpi=300, fmt_string='o'):
    """
    Creates a simple scatter plot of magnitude with time
    :param catalogue:
        Earthquake catalogue as instance of :class:
        hmtk.seismicity.catalogue.Catalogue
    :param bool plot_error:
        Choose to plot error bars (True) or not (False)
    :param str fmt_string:
        Symbology of plot
    """
    plt.figure(figsize=DEFAULT_SIZE)
    dtime = catalogue.get_decimal_time()
    if len(catalogue.data['sigmaMagnitude']) == 0:
        print 'Magnitude Error is missing - neglecting error bars!'
        plot_error = False

    if plot_error:
        plt.errorbar(dtime,
                     catalogue.data['magnitude'],
                     xerr=None,
                     yerr=catalogue.data['sigmaMagnitude'],
                     fmt=fmt_string)
    else:
        plt.plot(dtime, catalogue.data['magnitude'], fmt_string)
    plt.xlabel('Year', fontsize='large')
    plt.ylabel('Magnitude', fontsize='large')
    plt.title('Magnitude-Time Plot', fontsize='large')
    plt.show()
    _save_image(filename, filetype, dpi)
    return