def plot(self, outf=None, dosave=True, savedir="Plot/", show=True):
     if outf is None:
         outf = self.outf
         # print outf
     oo = mlab.csv2rec(outf, delimiter=" ")
     # print oo
     plt.errorbar(oo["time"] % self.period, oo["magnitude"], oo["error"], fmt="b.")
     plt.plot(oo["time"] % self.period, oo["model"], "ro")
     plt.title(
         "#%i P=%f d (chisq/dof = %f) r1+r2=%f"
         % (self.dotastro_id, self.period, self.outrez["chisq"], self.outrez.get("r1") + self.outrez.get("r2"))
     )
     ylim = plt.ylim()
     # print ylim
     if ylim[0] < ylim[1]:
         plt.ylim(ylim[1], ylim[0])
     plt.draw()
     if show:
         plt.show()
     if dosave:
         if not os.path.isdir(savedir):
             os.mkdir(savedir)
         plt.savefig("%splot%i.png" % (savedir, self.dotastro_id))  # ,self.period))
         print("Saved", "%splot%i.png" % (savedir, self.dotastro_id))  # ,self.period)
     plt.clf()

def display_vecstascollector_summary(stcol):
    n = stcol.getStat("N[0]")
    print n
    mean = stcol.getMean()
    stddev = stcol.getStdDev()
    ucov = (1.0/n)*stcol.getXtX()
    cov = stcol.getCovariance()
    corr = stcol.getCorrelation()
    pxy_px_py = ucov-outer(mean,mean)
    
    plt.subplot(4,2,1)
    plt.errorbar(arange(len(mean)),mean,yerr=stddev)
    plt.title("activations mean and stddev")
    plt.subplot(4,2,2)
    plot_histogram(mean, "activations mean")

    plt.subplot(4,2,3)
    plot_offdiag_histogram(ucov, "uncentered covariances")
    plt.subplot(4,2,4)
    plot_diag_histogram(ucov, "uncentered variances")

    plt.subplot(4,2,5)
    plot_offdiag_histogram(cov, "covariances")
    plt.subplot(4,2,6)
    plot_diag_histogram(cov, "variances")

    plt.subplot(4,2,7)
    plot_offdiag_histogram(corr, "correlations")
    plt.subplot(4,2,8)
    plot_histogram(stddev, "stddevs")

    plt.show()

 def get_fitting_isochrone(self, isochrone_fitter, ax=None,
                           spectrum_ids=[0,1]):
     combined_stellar_params = self.get_combined_stellar_parameters(
         spectrum_ids=spectrum_ids)
     (age, metallicity), closest_isochrone = \
         isochrone_fitter.find_closest_isochrone(
         teff=combined_stellar_params['teff'],
         logg=combined_stellar_params['logg'],
         feh=combined_stellar_params['feh'])
     if ax is None:
         return closest_isochrone
     else:
         ax.plot(closest_isochrone['teff'], closest_isochrone['logg'])
         ax.errorbar([combined_stellar_params['teff']],
                     [combined_stellar_params['logg']],
                     xerr=[combined_stellar_params['teff_uncertainty']],
                     yerr=[combined_stellar_params['logg_uncertainty']],
                     label='Teff={0:.2f}+-{1:.2f}\n logg={2:.2f}+-{3:.2f}'
                     .format(combined_stellar_params['teff'],
                             combined_stellar_params['teff_uncertainty'],
                             combined_stellar_params['logg'],
                             combined_stellar_params['logg_uncertainty']))
         ax.set_title('Age = {0:.2g} Gyr [Fe/H] = {1:.2g}'.format(age, metallicity))
         ax.invert_xaxis()
         ax.invert_yaxis()

def plot_runs(runs):
    """ Plot population evolutions
    """
    ts = range(len(runs[0]))
    cmap = plt.get_cmap('viridis')
    for i, r in enumerate(runs):
        mean, var = zip(*r)
        bm, cm = zip(*mean)
        bv, cv = zip(*var)

        color = cmap(float(i)/len(runs))

        plt.errorbar(ts, bm, fmt='-', yerr=bv, c=color)
        plt.errorbar(ts, cm, fmt='--', yerr=cv, c=color)

    plt.title('population evolution overview')
    plt.xlabel('time')
    plt.ylabel('value')

    plt.ylim((0, 1))

    plt.plot(0, 0, '-', c='black', label='benefit value')
    plt.plot(0, 0, '--', c='black', label='cost value')
    plt.legend(loc='best')

    plt.savefig('result.pdf')
    plt.show()

 def sanity_2dCircularFit(self):
   import numpy as np
   import matplotlib.pylab as plt
   from PyAstronomy import funcFit as fuf
   
   # Get the circular model and assign
   # parameter values
   c = fuf.Circle2d()
   c["r"] = 1.0
   c["t0"] = 0.0
   c["per"] = 3.0
   
   # Evaluate the model at a number of
   # time stamps
   t = np.linspace(0.0, 10.0, 20)
   pos = c.evaluate(t)
   
   # Add some error to the "measurement"
   pos += np.reshape(np.random.normal(0.0, 0.2, pos.size), pos.shape)
   err = np.reshape(np.ones(pos.size), pos.shape) * 0.2
   
   # Define free parameters and fit the model
   c.thaw(["r", "t0", "per"])
   c.fit(t, pos, yerr=err)
   c.parameterSummary()
   
   # Evaluate the model at a larger number of
   # points for plotting
   tt = np.linspace(0.0, 10.0, 200)
   model = c.evaluate(tt)
   
   # Plot the result
   plt.errorbar(pos[::,0], pos[::,1], yerr=err[::,1], \
                xerr=err[::,0], fmt='bp')
   plt.plot(model[::,0], model[::,1], 'r--')

def nova_plot():

	erg2mev=624151.

	fig=plot.figure()
	yrange = [1e-6,2e-4]
	xrange = [1e-1,1e5]
	plot.fill_between([0.2,10e3],[yrange[1],yrange[1]],[yrange[0],yrange[0]],facecolor='yellow',interpolate=True,color='yellow',alpha=0.5)
	plot.annotate('AMEGO',xy=(3,9e-5),xycoords='data',fontsize=26,color='black')

	lat=ascii.read("data/NMon2012.LAT.dat",names=['energy','en_low','en_high','flux','flux_err','tmp'])
	plot.scatter(lat['energy'],lat['flux']*erg2mev,color='red')
	plot.errorbar(lat['energy'],lat['flux']*erg2mev,xerr=[lat['en_low'],lat['en_high']],yerr=lat['flux_err']*erg2mev,ecolor='red',capsize=0,fmt='none')
	latul=ascii.read("data/NMon2012.LAT.limits.dat",names=['energy','en_low','en_high','flux','tmp1','tmp2','tmp3','tmp4'])
	plot.errorbar(latul['energy'],latul['flux']*erg2mev,xerr=[latul['en_low'],latul['en_high']],yerr=0.5*latul['flux']*erg2mev,uplims=True,ecolor='red',capsize=0,fmt='none')
	plot.scatter(latul['energy'],latul['flux']*erg2mev,color='red')

	leptonic=ascii.read("data/sp-NMon12-IC-best-fit-1MeV-30GeV.txt",names=['energy','flux'],data_start=1)
	hadronic=ascii.read("data/sp-NMon12-pi0-and-secondaries.txt",names=['energy','flux1','flux2'],data_start=1)	

	plot.plot(leptonic['energy'],leptonic['flux']*erg2mev,'r--',color='black',lw=2,label='Leptonic')
	plot.plot(hadronic['energy'],hadronic['flux2']*erg2mev,color='black',lw=2,label='Hadronic+Secondary Leptons')

	plot.legend(loc='upper right',fontsize='small',frameon=False,framealpha=0.5)
	plot.xscale('log')
	plot.yscale('log')
	plot.ylim(yrange)
	plot.xlim(xrange)
	plot.xlabel(r'Energy (MeV)')
	plot.ylabel(r'Energy$^2 \times $ Flux (Energy) (erg cm$^{-2}$ s$^{-1}$)')
	plot.title('Nova V339 Del 2013')
	plot.savefig('Nova_SED.png', bbox_inches='tight')
	plot.savefig('Nova_SED.eps', bbox_inches='tight')
	plot.show()
	plot.close()

def plot_avg(measure_type, tofile=None):
    fig,ax = plt.subplots()
    for name in system_names:
        Y = []
        syst, fmt = systems[name]
        for op in op_types:
            if measure_type == 'latency':
                latencies = get_datapoints_latency(op)
                Y_op, _ = latencies[name]
            elif measure_type == 'throughput':
                throughputs = get_datapoints_throughput(op)
                Y_op = throughputs[name]
            Y.append(Y_op)
        Y = np.mean(Y, axis=0)
        plt.errorbar(X,Y,fmt=fmt,label=name)
    ax.set_xticks(X)
    ax.set_xlabel('Number of concurrent nodes.')
    ax.set_xlim([0,17])
    if measure_type == 'throughput':
        ax.set_ylabel('Average throughput in KOps per second.')
    elif measure_type == 'latency':
        ax.set_ylabel('Average latency in ms.')
    lgd = plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
                     ncol=4, mode="expand", borderaxespad=0.)
    plt.grid(True)
    if tofile is not None:
        plt.savefig(tofile, bbox_extra_artists=(lgd,), bbox_inches='tight')
    else:
        plt.show()

def plotRocCurves(file_legend):
	pylab.clf()
	pylab.figure(1)
	pylab.xlabel('1 - Specificity', fontsize=12)
	pylab.ylabel('Sensitivity', fontsize=12)
	pylab.title("Need for Referral")
	pylab.grid(True, which='both')
	pylab.xticks([i/10.0 for i in range(1,11)])
	pylab.yticks([i/10.0 for i in range(0,11)])
	pylab.tick_params(axis="both", labelsize=15)

	for file, legend in file_legend:
		points = open(file,"rb").readlines()
		x = [float(p.split()[0]) for p in points]
		y = [float(p.split()[1]) for p in points]
		dev = [float(p.split()[2]) for p in points]
		x = [0.0] + x
		y = [0.0] + y
		dev = [0.0] + dev
	
		auc = np.trapz(y, x) * 100
		aucDev = np.trapz(dev, x) * 100

		pylab.grid()
		pylab.errorbar(x, y, yerr = dev, fmt='-')
		pylab.plot(x, y, '-', linewidth = 1.5, label = legend + u" (AUC = {0:0.1f}% \xb1 {1:0.1f}%)".format(auc,aucDev))

	pylab.legend(loc = 4, borderaxespad=0.4, prop={'size':12})
	pylab.savefig("referral/referral-curves.pdf", format='pdf')

    def do_plot_obs(self):
        """ Plot the observed radial velocities as a function of time.
        Data from each file are color coded and labeled.
        """
        # import pyqtgraph as pg

        colors = 'bgrcmykw' # lets hope for less than 9 data-sets
        t, rv, err = self.time, self.vrad, self.error # temporaries
        
        plt.figure()
        # p = pg.plot()
        # plot each files' values
        for i, (fname, [n, nout]) in enumerate(sorted(self.provenance.iteritems())):
            m = n-nout # how many values are there after restriction
            
            # e = pg.ErrorBarItem(x=t[:m], y=rv[:m], \
            #                     height=err[:m], beam=0.5,\
            #                     pen=pg.mkPen(None))
                                # pen={'color': 0.8, 'width': 2})
            # p.addItem(e)
            # p.plot(t[:m], rv[:m], symbol='o')
            plt.errorbar(t[:m], rv[:m], yerr=err[:m], \
                         fmt='o'+colors[i], label=fname)
            t, rv, err = t[m:], rv[m:], err[m:]
        
        plt.xlabel('Time [days]')
        plt.ylabel('RV [km/s]')
        plt.legend()
        plt.tight_layout()
        plt.show()

def plot(scores, scores2=None):
    import matplotlib.pylab as pl
    from matplotlib.ticker import FuncFormatter

    def percentages(x, pos=0):
        return "%2.2f%%" % (100 * x)

    ax1 = pl.subplot(211)
    pl.errorbar(scores2[:, 1], scores2[:, 2], yerr=scores2[:, 5], c="k", marker="o")
    # if scores2 is not None:
    #    pl.errorbar(scores2[:, 1] + 0.02, scores2[:, 2], yerr=scores2[:, 5],
    #            c='0.5', marker='s')
    pl.ylabel("Singular acc.")
    ax1.yaxis.set_major_formatter(FuncFormatter(percentages))

    pl.xlabel("Proportion of training set used")
    ax2 = pl.subplot(212, sharex=ax1)
    pl.errorbar(scores[:, 1], scores[:, 3], yerr=scores[:, 6], c="k", marker="o")
    if scores2 is not None:
        pl.errorbar(scores2[:, 1], scores2[:, 3], yerr=scores2[:, 6], c="k", marker="s")

    ax2.yaxis.set_major_formatter(FuncFormatter(percentages))

    # ax3 = pl.subplot(313, sharex=ax2)
    pl.errorbar(scores[:, 1] + 0.02, scores[:, 4], yerr=scores[:, 7], c="0.5", marker="o")
    if scores2 is not None:
        pl.errorbar(scores2[:, 1] + 0.02, scores2[:, 4], yerr=scores2[:, 7], c="0.5", marker="s")
    pl.ylabel("Plural and combined acc.")
    # ax3.yaxis.set_major_formatter(FuncFormatter(percentages))
    # pl.setp(ax3.get_xticklabels(), visible=False)

    # pl.show()

    for ext in ("pdf", "svg", "png"):
        pl.savefig("train_size-i.%s" % ext)

        def parse_error(d):
            if not np.isfinite(d.fp_local):
                return
            f = ais_code.fp_to_ne(d.fp_local)
            f0 = ais_code.fp_to_ne(d.fp_local + d.fp_local_error)
            f1 = ais_code.fp_to_ne(d.fp_local - d.fp_local_error)
            if errors:
                plt.plot((d.time, d.time),(f0,f1),
                    color='lightgrey', linestyle='solid',
                    marker='None', zorder=-1000,**kwargs)
            plt.plot(d.time, f, fmt, ms=self.marker_size, zorder=1000, **kwargs)

            if full_marsis and hasattr(d, 'maximum_fp_local'):
                plt.plot(d.time, ais_code.fp_to_ne(d.maximum_fp_local),
                    'b.', ms=self.marker_size, zorder=900, **kwargs)

            if full_marsis:
                if np.isfinite(d.morphology_fp_local):
                    v, e = ais_code.fp_to_ne(d.morphology_fp_local,
                                                d.morphology_fp_local_error)
                    plt.errorbar(float(d.time), v, yerr=e,
                            marker='x', ms=1.3, color='purple',
                            zorder=1e90, capsize=0., ecolor='plum')

                if np.isfinite(d.integrated_fp_local):
                    v, e = ais_code.fp_to_ne(d.integrated_fp_local,
                                                d.integrated_fp_local_error)
                    plt.errorbar(float(d.time), v, yerr=e,
                            marker='x', ms=1.3, color='blue',
                            zorder=1e99, capsize=0., ecolor='cyan')

def errorbar_plot(data, x_spec, y_spec, fname):
    """ Dynamically create errorbar plot
    """
    x_label, x_func = x_spec
    y_label, y_func = y_spec

    # compute data
    points = collections.defaultdict(list)
    for syst, mat, _ in data:
        x_value = x_func(syst, mat)
        y_value = y_func(syst, mat)

        if x_value is None or y_value is None: continue
        points[x_value].append(y_value)

    # plot figure
    densities = []
    averages = []
    errbars = []
    for dens, avgs in points.items():
        densities.append(dens)
        averages.append(np.mean(avgs))
        errbars.append(np.std(avgs))

    plt.errorbar(
        densities, averages, yerr=errbars,
        fmt='o', clip_on=False)

    plt.title('')
    plt.xlabel(x_label)
    plt.ylabel(y_label)

    plt.tight_layout()
    save_figure('images/%s' % fname, bbox_inches='tight')
    plt.close()

def compareMassRatioVsRedshift(his, hiserr, mine, myerr, hisfield, hisid, myfield, myid, z):   
    import matplotlib.pylab as plt
         
    a = []
    b = []
    c = []
    d = []
    e = []
    z_list = []
    
    for i in range(0, np.shape(mine)[0]):
        if mine[i] != 0:
            #his_id_idx = np.where(hisid == myid[i])[0]
            pos_id_idx = np.where(hisid == myid[i])[0]
            if len(pos_id_idx) == 1:
                his_id_idx = int(pos_id_idx)
            else:
                pos_field_idx = np.where(hisfield == myfield[i])[0]
                for ii in range(0, len(pos_field_idx)):
                    for iii in range(0, len(pos_id_idx)):
                        if pos_field_idx[ii] == pos_id_idx[iii]:
                            his_id_idx = int(pos_id_idx[iii])
    
            #print str(myid[i]), str(myfield[i]), str(hisid[his_id_idx]), str(hisfield[his_id_idx])
            assert myfield[i] == hisfield[his_id_idx]
            assert myid[i] == hisid[his_id_idx]
            if his[his_id_idx] != -999:
                a.append(his[his_id_idx])
                b.append(hiserr[his_id_idx])
                c.append(mine[i])#-np.log10(((1.0+z[i]))))
                d.append(myerr[i,0])#-np.log10(((1.0+z[i]))))
                e.append(myerr[i,1])#-np.log10(1.0+z[i]))
                z_list.append(z[i])
                
    # since errors are just percentile, need difference from median
    d = np.asarray(c)-np.asarray(d)
    e = np.asarray(e)-np.asarray(c)
    
    c_a = 10**np.array(c)
    a_a = 10**np.array(a)
    
    ratio = c_a/a_a
    
    plt.errorbar(z_list, ratio, fmt='o', 
                 color='b', capsize=0, alpha=0.50)
                 
    # plot the y=x line
    x = np.linspace(np.min(z_list),
                    np.max(z_list),
                    10)
    plt.plot(x, np.ones(len(x)), 'k--')
                 
    plt.yscale('log')
    
    plt.xlabel("Redshift")
    plt.ylabel("MCSED  / Alex's  [note: ratio of actual masses]")
    plt.title("Redshift vs Mass Ratio (MCSED/Alex)")
    #plt.legend(['With Neb. Emis.'],loc=0)#, 'W/o Neb. Emis'], loc=0)
        
    plt.show()

def main():
    parser = OptionParser(description='Fitting to a noisy data generated by a known function')
    parser.add_option("--npoints", type="int",   help="number of data points") 
    parser.add_option("--low",     type="float", help="smallest data point") 
    parser.add_option("--high",    type="float", help="highest data point") 
    parser.add_option("--sigma",   type="float", help="std of noise") 
    (options, args) = parser.parse_args() 

    pl.figure(1,(7,6))
    ax = pl.subplot(1,1,1)

    pl.connect('key_press_event',kevent.press)
    
    
    sigma = options.sigma    
    Ls   = np.append(np.linspace(options.low,options.high,options.npoints),46)
    nLs  = np.linspace(min(Ls),max(Ls),100)
    Mis  = HalfLog(Ls,.5,0.5)
    errs = np.random.normal(0,sigma, len(Mis))
    Mis  = Mis+errs
    pl.errorbar(Ls,Mis,errs,ls='',marker='s',color='b')
    print sigma/Mis 

    coeff, var_matrix = curve_fit(FreeLog,Ls,Mis,(1.0,1.0,1.0))
    err = np.sqrt(np.diagonal(var_matrix))
    dof     = len(Ls) - len(coeff)
    chisq   = sum(((Mis-FreeLog(Ls,coeff[0],coeff[1],coeff[2]))/sigma)**2)
    cdf     = special.chdtrc(dof,chisq)
    print 'Free:  a = %0.2f(%0.2f); b = %0.2f(%0.2f); c = %0.2f(%0.2f); p-value = %0.2f ' %(coeff[0],err[0],coeff[1],err[1],coeff[2],err[2],cdf)
    pl.plot(nLs,FreeLog(nLs,coeff[0],coeff[1],coeff[2]),label='Free',color='y')

    coeff, var_matrix = curve_fit(ZeroLog,Ls,Mis,(1.0,1.0))
    err = np.sqrt(np.diagonal(var_matrix))
    dof     = len(Ls) - len(coeff)
    chisq   = sum(((Mis-ZeroLog(Ls,coeff[0],coeff[1]))/sigma)**2)
    cdf     = special.chdtrc(dof,chisq)
    print 'Zero:  a = %0.2f(%0.2f);                 c = %0.2f(%0.2f); p-value = %0.2f' %(coeff[0],err[0],coeff[1],err[1],cdf)
    pl.plot(nLs,ZeroLog(nLs,coeff[0],coeff[1]),label='Zero',color='g')
    pl.tight_layout()

    coeff, var_matrix = curve_fit(HalfLog,Ls,Mis,(1.0,1.0))
    err = np.sqrt(np.diagonal(var_matrix))
    dof     = len(Ls) - len(coeff)
    chisq   = sum(((Mis-HalfLog(Ls,coeff[0],coeff[1]))/sigma)**2)
    cdf     = special.chdtrc(dof,chisq)
    print 'Half:  a = %0.2f(%0.2f);                 c = %0.2f(%0.2f); p-value = %0.2f' %(coeff[0],err[0],coeff[1],err[1],cdf)
    pl.plot(nLs,HalfLog(nLs,coeff[0],coeff[1]),label='Half',color='b')
    pl.tight_layout()

    coeff, var_matrix = curve_fit(OneLog,Ls,Mis,(1.0,1.0))
    err = np.sqrt(np.diagonal(var_matrix))
    dof     = len(Ls) - len(coeff)
    chisq   = sum(((Mis-OneLog(Ls,coeff[0],coeff[1]))/sigma)**2)
    cdf     = special.chdtrc(dof,chisq)
    print 'Unity: a = %0.2f(%0.2f);                 c = %0.2f(%0.2f); p-value = %0.2f' %(coeff[0],err[0],coeff[1],err[1],cdf)
    pl.plot(nLs,OneLog(nLs,coeff[0],coeff[1]),label='Unity',color='r')
    pl.tight_layout()

    pl.legend()
    pl.show()

def extractEffectiveMass(df,cut_low=0,cut_right=None,n_cuts=5,makePlot=False):
    if cut_right==None:
        cut_right=max(df["t"])
    window=(cut_right-cut_low)/n_cuts
    
    slope=np.zeros(n_cuts)
    intercept=np.zeros(n_cuts)
    slopeError=np.zeros(n_cuts)
    interceptError=np.zeros(n_cuts)
    
    for i in range(0,n_cuts):
        # select the right intervals for the linear fit
        cut_low_current=cut_low + i*window
        cut_high_current=cut_low + (i+1)*window
        df1=df[(df["t"]>cut_low_current) & (df["t"]<cut_high_current) ]
        if len(df1["t"]) <= 3:
            raise not_enough_data()
        
        params,covs=curve_fit(linear,df1["t"],df1["W"],sigma=df1["deltaW"],maxfev=100000)
        slope[i]=params[0]
        intercept[i]=params[1]
        slopeError[i]=sqrt(covs[0][0])
        interceptError[i]=sqrt(covs[1][1])
        
        if makePlot==True:
            up=(slope[i]+slopeError[i])+(intercept[i]+interceptError[i])/df1["t"]
            down=(slope[i]-slopeError[i])+(intercept[i]-interceptError[i])/df1["t"]
            
            plt.fill_between(df1["t"],up,down,alpha=0.4)
            plt.errorbar(np.array(df1["t"]),np.array(df1["W"])/np.array(df1["t"]),np.array(df1["deltaW"])/np.array(df1["t"]),fmt="or")
    return np.array([slope.mean(),sqrt( slopeError.mean()**2 + slope.var())])

 def sanity_example_binningx0dt_example2(self):
   """
     Checking `binningx0dt` example 2.
   """
   import numpy as np
   import matplotlib.pylab as plt
   from PyAstronomy.pyasl import binningx0dt
   
   # Generate some data
   x = np.arange(-100,999)
   # Create some holes in the data
   x = np.delete(x, range(340,490))
   x = np.delete(x, range(670,685))
   x = np.delete(x, range(771,779))
   y = np.sin(x/100.)
   y += np.random.normal(0,0.1,len(x))
   
   
   # Bin using bin width of 27 and starting at minimum x-value.
   # Use beginning of bin as starting value.
   r1, dt1 = binningx0dt(x, y, dt=27, x0=min(x), useBinCenter=True)
   
   # As previously, but use the mean x-value in the bins to produce the
   # rebinned time axis.
   r2, dt2 = binningx0dt(x, y, dt=27, x0=min(x), useMeanX=True)
   
   print "Median shift between the time axes: ", np.median(r1[::,0] - r2[::,0])
   print " -> Time bins are not aligned due to 'forced' positioning of"
   print "    the first axis."
   
   # Plot the output
   plt.plot(x,y, 'b.-')
   plt.errorbar(r1[::,0], r1[::,1], yerr=r1[::,2], fmt='kp--')
   plt.errorbar(r2[::,0], r2[::,1], yerr=r2[::,2], fmt='rp--')

def method_errorbar(data,xlabels, line_color=default_color,
                   med_color=None, legend=None, y_offset=0.0,alpha=0.05):
    if not med_color:
        med_color=line_color
    ax.grid(axis='x', color='0.9', linestyle='-', linewidth=0.2)
    ax.set_axisbelow(True)
    n,m=data.shape
    medians = [percentile(data[:,i],50) for i in range(m)]
    xerr = [[ medians[i]-percentile(data[:,i],100*(alpha/2.)),
              percentile(data[:,i],100*(1-alpha/2.))-medians[i] ]
             for i in range(m)]
    xerr=np.array(xerr).transpose()
    y_marks = np.array(range(len(xlabels)))-y_offset
    plt.errorbar(y=y_marks,
                 x=medians,xerr=xerr,fmt='|',capsize=0,color=line_color,
                 ecolor=line_color,elinewidth=0.3,markersize=2)
    plt.xlabel('% cases used', fontsize=8)
    ax.tick_params(axis='x', which='both', labelsize=8)
    ax.set_yticks(np.array(range(len(xlabels))))
    ax.set_yticklabels(xlabels,fontsize=6)
    plt.ylim((min(y_marks)-0.5,max(y_marks)+0.5))
    spines_to_remove = ['top', 'right','left']
    for spine in spines_to_remove:
        ax.spines[spine].set_visible(False)
    ppl.utils.remove_chartjunk(ax, ['top', 'right', 'bottom'], show_ticks=False)
    if legend:
        rect = legend.get_frame()
        rect.set_facecolor(light_grey)
        rect.set_linewidth(0.0)

def plot_source_radius(cat):
    mw = []
    mw_std = []
    source_radius = []
    source_radius_std = []

    plt.figure(figsize=(10, 4.5))

    # Read the source radius.
    for event in cat:
        mag = event.magnitudes[1]
        if len(mag.comments) != 2:
            continue
        mw.append(mag.mag)
        mw_std.append(mag.mag_errors.uncertainty)
        sr, std = mag.comments[1].text.split(";")
        _, sr = sr.split("=")
        _, std = std.split("=")
        sr = float(sr[:-1])
        std = float(std)
        source_radius.append(sr)
        source_radius_std.append(std)
    plt.errorbar(mw, source_radius, yerr=source_radius_std,
        fmt="o", linestyle="None")
    plt.xlabel("Mw", fontsize="x-large")
    plt.ylabel("Source Radius [m]", fontsize="x-large")
    plt.grid()
    plt.savefig("/Users/lion/Desktop/SourceRadius.pdf")

def plot_table(*positional_parameters,errorbars=None,xrange=None,yrange=None,
               title="",xtitle="",ytitle="",show=1,legend=None,color=None):

    n_arguments = len(positional_parameters)
    if n_arguments == 0:
        return

    fig = plt.figure()

    if n_arguments == 1:
        y = positional_parameters[0]
        x = np.arange(y.size)
        plt.plot(x,y,label=legend)
    elif n_arguments == 2:
        x = positional_parameters[0]
        y = positional_parameters[1]

        if len(y.shape) == 1:
            y = np.reshape(y,(1,y.size))
            if isinstance(legend,str):
                legend = [legend]
            if isinstance(color,str):
                color = [color]

        for i in range(y.shape[0]):
            if legend is None:
                ilegend = None
            else:
                ilegend = legend[i]

            if color is None:
                icolor = None
            else:
                icolor = color[i]

            if errorbars is None:
                plt.plot(x,y[i],label=ilegend,color=icolor)
            else:
                plt.errorbar(x,y[i],yerr=errorbars[i],label=ilegend,color=icolor)
    else:
        raise Exception("Incorrect number of arguments")


    plt.xlim( xrange )
    plt.ylim( yrange )

    if legend is not None:
        ax = plt.subplot(111)
        ax.legend(bbox_to_anchor=(1.1, 1.05))

    plt.title(title)
    plt.xlabel(xtitle)
    plt.ylabel(ytitle)


    if show:
        plt.show()

    return fig

def main():
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument('losses', nargs='+', type=str)
    parser.add_argument('-o', '--output', default='log-loss', type=str)
    #parser.add_argument('--dataset', default='Test', type=str)
    parser.add_argument('--captions', nargs='+', type=str)

    args = parser.parse_args()

    rs = np.random.RandomState(0)

    vz = VzLog(args.output)

    plt.figure()
    for i, loss_fn in enumerate(args.losses):
        print('file', loss_fn)
        data = dd.io.load(loss_fn)
        iter = data['iterations'][0]
        losses = data['losses'].mean(0)
        losses_std = data['losses'].std(0)

        if args.captions:
            caption = args.captions[i]
        else:
            caption = loss_fn

        caption = '{} ({:.3f})'.format(caption, losses[-1])

        plt.errorbar(iter, losses, yerr=losses_std, label=caption)

    plt.legend()
    plt.ylabel('Loss')
    plt.xlabel('Iteration')
    vz.savefig()
    plt.close()

    plt.figure()
    for i, loss_fn in enumerate(args.losses):
        print('file', loss_fn)
        data = dd.io.load(loss_fn)
        iter = data['iterations'][0]
        rates = 100*(1-data['rates'].mean(0))
        rates_std = 100*data['rates'].std(0)

        if args.captions:
            caption = args.captions[i]
        else:
            caption = loss_fn

        caption = '{} ({:.2f}%)'.format(caption, rates[-1])

        plt.errorbar(iter, rates, yerr=rates_std, label=caption)

    plt.legend()
    plt.ylabel('Error rate (%)')
    plt.xlabel('Iteration')
    vz.savefig()
    plt.close()