def test_power_spectra(r0, N, delta, L0, l0):
    N*= 10
    phase_screen = atmosphere.ft_phase_screen(r0, N, delta, L0, l0)
    phase_screen = phase_screen[:N/10, :N/10]
    power_spec_2d = numpy.fft.fft2(phase_screen, s=(N*2, N*2))

    plt.figure()
    plt.imshow(numpy.abs(numpy.fft.fftshift(power_spec_2d)), interpolation='nearest')

    power_spec = circle.aziAvg(numpy.abs(numpy.fft.fftshift(power_spec_2d)))
    power_spec /= power_spec.sum()

    freqs = numpy.fft.fftfreq(power_spec_2d.shape[0], delta)

    # Theoretical Model of Power Spectrum

    print freqs

    plt.figure()
    plt.plot(freqs[:freqs.size/2], power_spec)
    plt.xscale('log')
    plt.yscale('log')
    plt.show()

    return None

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 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_ccdf(values, xscale, yscale):
    pylab.yscale(yscale)
    cdf = Cdf.MakeCdfFromList(values)
    values, prob = cdf.Render()
    pylab.xscale(xscale)
    compProb = [1 - e for e in prob] 
    pylab.plot(values, compProb)

def smooth_color_prior(size=64, sigma=5, do_plot=False):

    prior_prob = np.load(os.path.join(data_dir, "CelebA_%s_prior_prob.npy" % size))
    # add an epsilon to prior prob to avoid 0 vakues and possible NaN
    prior_prob += 1E-3 * np.min(prior_prob)
    # renormalize
    prior_prob = prior_prob / (1.0 * np.sum(prior_prob))

    # Smooth with gaussian
    f = interp1d(np.arange(prior_prob.shape[0]),prior_prob)
    xx = np.linspace(0,prior_prob.shape[0] - 1, 1000)
    yy = f(xx)
    window = gaussian(2000, sigma)  # 2000 pts in the window, sigma=5
    smoothed = convolve(yy, window / window.sum(), mode='same')
    fout = interp1d(xx,smoothed)
    prior_prob_smoothed = np.array([fout(i) for i in range(prior_prob.shape[0])])
    prior_prob_smoothed = prior_prob_smoothed / np.sum(prior_prob_smoothed)

    # Save
    file_name = os.path.join(data_dir, "CelebA_%s_prior_prob_smoothed.npy" % size)
    np.save(file_name, prior_prob_smoothed)

    if do_plot:
        plt.plot(prior_prob)
        plt.plot(prior_prob_smoothed, "g--")
        plt.plot(xx, smoothed, "r-")
        plt.yscale("log")
        plt.show()

def plotFeaturePDF(ift, pft, outbase, fmin=0.0, fmax=1.0, fstep=0.01):
    """
    Plot a comparison between the input feature distribution and the 
    feature distribution of the predicted halos
    """
    plt.clf()
    nfbins = ( fmax - fmin ) / fstep
    fbins = np.logspace( fmin, fmax, nfbins )
    fcen = ( fbins[:-1] + fbins[1:] ) / 2

    plt.xscale( 'log', nonposx='clip' )
    plt.yscale( 'log', nonposy='clip' )
    
    ic, e, p = plt.hist( ift, fbins, label='Original Halos', alpha=0.5, normed=True )
    pc, e, p = plt.hist( pft, fbins, label='Added Halos', alpha=0.5, normed=True )

    plt.legend()
    plt.xlabel( r'$\delta$' )
    plt.savefig( outbase+'_fpdf.png' )

    fdtype = np.dtype( [ ('fcen', float), ('ifcounts', float), ('pfcounts', float) ] )
    fd = np.ndarray( len(fcen), dtype = fdtype )
    fd[ 'mcen' ] = fcen
    fd[ 'imcounts' ] = ic
    fd[ 'pmcounts' ] = pc

    fitsio.write( outbase+'_fpdf.fit', fd )

def plotMassFunction(im, pm, outbase, mmin=9, mmax=13, mstep=0.05):
    """
    Make a comparison plot between the input mass function and the 
    predicted projected correlation function
    """
    plt.clf()

    nmbins = ( mmax - mmin ) / mstep
    mbins = np.logspace( mmin, mmax, nmbins )
    mcen = ( mbins[:-1] + mbins[1:] ) /2
    
    plt.xscale( 'log', nonposx = 'clip' )
    plt.yscale( 'log', nonposy = 'clip' )
    
    ic, e, p = plt.hist( im, mbins, label='Original Halos', alpha=0.5, normed = True)
    pc, e, p = plt.hist( pm, mbins, label='Added Halos', alpha=0.5, normed = True)
    
    plt.legend()
    plt.xlabel( r'$M_{vir}$' )
    plt.ylabel( r'$\frac{dN}{dM}$' )
    #plt.tight_layout()
    plt.savefig( outbase+'_mfcn.png' )
    
    mdtype = np.dtype( [ ('mcen', float), ('imcounts', float), ('pmcounts', float) ] )
    mf = np.ndarray( len(mcen), dtype = mdtype )
    mf[ 'mcen' ] = mcen
    mf[ 'imcounts' ] = ic
    mf[ 'pmcounts' ] = pc

    fitsio.write( outbase+'_mfcn.fit', mf )

def get_OH(OIII4363,OIII4959,OIII5007,Hb):
	Te = np.arange(5000,20000,1)
	t3 = Te/1e4
	ne = 100	# cm^-3
	x = 1e-4*ne*t3**(-0.5)
	C_T = (8.44-1.09*t3+0.5*t3**2.-0.08*t3**3.)*(1.+0.0004*x)/(1.+0.044*x)
	log_OIII_ratio = 1.432/t3+np.log10(C_T)
	log_OIII_ratio_obs = np.log10((OIII4959+OIII5007)/OIII4363)

	Te_obs = []
	plt.clf()
	plt.plot(Te,log_OIII_ratio,color='black',marker='.',linestyle='none')
	plt.yscale('log')
	for i in range(len(OIII4363)):
		plt.axhline(log_OIII_ratio_obs[i],linestyle='--')
		d_ratio = abs(log_OIII_ratio_obs[i]-log_OIII_ratio)
		min_d_ratio = min(d_ratio)
		min_sub = list(d_ratio).index(min_d_ratio)
		Te_obs.append(Te[min_sub])
	plt.xlim(10000,20000)
	plt.ylim(1.,3)
	plt.xlabel('Te')
	plt.ylabel('(OIII4959+5007)/OIII4363')

	Te_obs = np.array(Te_obs)
	t3_obs = Te_obs/1e4
	logOIIIH = np.log10((OIII4959+OIII5007)/Hb)+6.200+1.251+1.251/t3_obs - \
				5*np.log10(t3_obs)-0.014*t3_obs
	t2_obs = -0.577+t3*(2.065-0.498*t3)
	logOIIH = np.log10(OII3727/Hb)+5.961+1.676/t2_obs-0.4*np.log10(t2_obs) - \
				0.034*t2_obs+np.log10(1+1.35*x)
	OH = 10**(logOIIIH-12.)+10**(logOIIIH-12.)
	logOH = 12 + np.log10(OH)

	return Te_obs,logOIIH,logOIIIH,logOH

 def test_simple_gen(self):
     self_con = .8
     other_con = 0.05
     g = self.gen.gen_stoch_blockmodel(min_degree=1, blocks=5, self_con=self_con, other_con=other_con,
                                       powerlaw_exp=2.1, degree_seq='powerlaw', num_nodes=1000, num_links=3000)
     deg_hist = vertex_hist(g, 'total')
     res = fit_powerlaw.Fit(g.degree_property_map('total').a, discrete=True)
     print 'powerlaw alpha:', res.power_law.alpha
     print 'powerlaw xmin:', res.power_law.xmin
     if len(deg_hist[0]) != len(deg_hist[1]):
         deg_hist[1] = deg_hist[1][:len(deg_hist[0])]
     print 'plot degree dist'
     plt.plot(deg_hist[1], deg_hist[0])
     plt.xscale('log')
     plt.xlabel('degree')
     plt.ylabel('#nodes')
     plt.yscale('log')
     plt.savefig('deg_dist_test.png')
     plt.close('all')
     print 'plot graph'
     pos = sfdp_layout(g, groups=g.vp['com'], mu=3)
     graph_draw(g, pos=pos, output='graph.png', output_size=(800, 800),
                vertex_size=prop_to_size(g.degree_property_map('total'), mi=2, ma=30), vertex_color=[0., 0., 0., 1.],
                vertex_fill_color=g.vp['com'],
                bg_color=[1., 1., 1., 1.])
     plt.close('all')
     print 'init:', self_con / (self_con + other_con), other_con / (self_con + other_con)
     print 'real:', gt_tools.get_graph_com_connectivity(g, 'com')

def compute_color_prior(size=64, do_plot=False):

    # Load the gamut points location
    q_ab = np.load(os.path.join(data_dir, "pts_in_hull.npy"))

    if do_plot:
        plt.figure(figsize=(15, 15))
        gs = gridspec.GridSpec(1, 1)
        ax = plt.subplot(gs[0])
        for i in range(q_ab.shape[0]):
            ax.scatter(q_ab[:, 0], q_ab[:, 1])
            ax.annotate(str(i), (q_ab[i, 0], q_ab[i, 1]), fontsize=6)
            ax.set_xlim([-110,110])
            ax.set_ylim([-110,110])

    with h5py.File(os.path.join(data_dir, "CelebA_%s_data.h5" % size), "a") as hf:
        # Compute the color prior over a subset of the training set
        # Otherwise it is quite long
        X_ab = hf["training_lab_data"][:100000][:, 1:, :, :]
        npts, c, h, w = X_ab.shape
        X_a = np.ravel(X_ab[:, 0, :, :])
        X_b = np.ravel(X_ab[:, 1, :, :])
        X_ab = np.vstack((X_a, X_b)).T

        if do_plot:
            plt.hist2d(X_ab[:, 0], X_ab[:, 1], bins=100, norm=LogNorm())
            plt.xlim([-110, 110])
            plt.ylim([-110, 110])
            plt.colorbar()
            plt.show()
            plt.clf()
            plt.close()

        # Create nearest neighbord instance with index = q_ab
        NN = 1
        nearest = nn.NearestNeighbors(n_neighbors=NN, algorithm='ball_tree').fit(q_ab)
        # Find index of nearest neighbor for X_ab
        dists, ind = nearest.kneighbors(X_ab)

        # We now count the number of occurrences of each color
        ind = np.ravel(ind)
        counts = np.bincount(ind)
        idxs = np.nonzero(counts)[0]
        prior_prob = np.zeros((q_ab.shape[0]))
        for i in range(q_ab.shape[0]):
            prior_prob[idxs] = counts[idxs]

        # We turn this into a color probability
        prior_prob = prior_prob / (1.0 * np.sum(prior_prob))

        # Save
        np.save(os.path.join(data_dir, "CelebA_%s_prior_prob.npy" % size), prior_prob)

        if do_plot:
            plt.hist(prior_prob, bins=100)
            plt.yscale("log")
            plt.show()

def plot_change(alldata,thisdata,name):
	pl.style.use('ggplot')
	f, ax = pl.subplots()
	
	ax.loglog(alldata,'.',label='All mirNA collected in DB')
	ax.loglog(alldata.index(thisdata),thisdata,'o',ms=9,label=name)	
	ax.legend(loc=2)
	pl.ylabel('Expression change')
	pl.yscale('log')
	pl.ylim(alldata[2],alldata[-2])
	pl.yticks(np.logspace(np.log10(alldata[2]),np.log10(alldata[-2]),10),[i.round(2) for i in np.logspace(np.log10(alldata[2]),np.log10(alldata[-2]),10)])

	canvas = FigureCanvas(f)
	response = HttpResponse(content_type='image/png')
	canvas.print_png(response)
	return response

    def plot_noise(self, apsize):

        logging.info('Plotting noise budget')

        p = self.p.aperture[apsize]

        plt.scatter(p.medmag, p.rms)
        plt.plot(p.mags_model, p.sigmastar_model, label='Star noise')
        plt.plot(p.mags_model, p.sigmasky_model, label='Sky noise')
        plt.plot(p.mags_model, p.sigmaread_model, label='Readout noise')
        plt.plot(p.mags_model, p.sigscint_model, label='Scint. noise')
        plt.yscale('log')
        plt.plot(p.mags_model, p.sigtot_model, label='Tot. noise')
        plt.legend(loc=4)

        plt.savefig(os.path.join(self.plot_out, 'noise_budget.pdf'))
        plt.show()

def plotPSD(lcInt, shortExp,**kwargs):
    '''
    plot power spectral density of lc
    return frequencies and powers from periodogram
    '''
    freq = 1.0/shortExp
    f, p = periodogram(lcInt,fs = 1./shortExp)

    plt.plot(f,p/np.max(p),**kwargs)
    plt.xlabel(r"Frequency (Hz)",fontsize=14)
    plt.xscale('log')
    plt.ylabel(r"Normalized PSD",fontsize=14)
    plt.yscale('log')
    plt.title(r"Lightcurve Power Spectrum",fontsize=14)
    plt.show()

    return f,p

def show_reverseattn(trange_gaincal, trange_other=None, first_attn_base=5, second_attn_base=3, corrected_attns=False):
        # This routine returns the same things as get_reverseattn, but it also graphs the three different "tsys"
        #  variables. It interactively asks the user which antenna, which polarization, and which frequency they 
        #  wish to view. Once it creates a graph, it will ask the user whether they wish to see more graphs or not.
        #  The first parameter it takes should be a GAINCALTEST trange, which may be located with find_gaincal, 
        #  The second parameter it takes should be a FEATTNTEST trange, which may be located with find_gaincal as
        #  well. It may or may not take a third parameter. Any form of file that was recorded by the system from 
        #  from the antennas may be inserted here, whether a flare or just quiet sun data. 
        #
        #  PLEASE NOTE: ANY trange with data may be used as "trange_gaincal", use a trange from a GAINCALTEST and the other 
        #  file as "trange_other" if you want the noise to be calibrated from the GAINCALTEST file, which will most likely
        #  be more recent than the FEATTNTEST file it would otherwise take the noise from. 
        tsys_attn_noise_corrected, tsys_noise_corrected, tsys = get_reverseattn(trange_gaincal, trange_other, first_attn_base, second_attn_base, corrected_attns)
        plotting_ = True
        while plotting_ == True:
            print ' '
            antenna = input('Which antenna would you like to see? You can choose from 1 to 8 : ')
            print ' '
            polarization = input('Which polarization would you like to see? Enter 0 for x, and 1 for y: ') 
            print ' '
            channel = input('Which channel would you like to see? You can choose from channel 0 to ' +str(tsys_attn_noise_corrected.shape[2]) + ' : ')

            antenna = int(antenna)-1
            polarization = int(polarization)
            channel = int(channel)

            plt.figure()
            tsys_ = plt.plot(tsys[antenna, polarization, channel, :], label='tsys')
            tsys_noise_corrected_ = plt.plot(tsys_noise_corrected[antenna, polarization, channel, :], label='noise corrected tsys')
            tsys_attn_noise_corrected_ = plt.plot(tsys_attn_noise_corrected[antenna, polarization, channel, :], label='noise and attn corrected tsys')
            plt.yscale('log')
            fontP = FontProperties()
            fontP.set_size('small')
            plt.legend(prop = fontP)
            plt.show()

            print ' '
            yesorno = input('Would you like to see another plot? Please answer 0 for yes or 1 for no. ')
            if yesorno == 0:
                pass
            else:
                print ' '
                print 'See you next time!'
                plotting_ = False
                return tsys_attn_noise_corrected, tsys_noise_corrected, tsys

def compute_prior_factor(size=64, gamma=0.5, alpha=1, do_plot=False):

    file_name = os.path.join(data_dir, "CelebA_%s_prior_prob_smoothed.npy" % size)
    prior_prob_smoothed = np.load(file_name)

    u = np.ones_like(prior_prob_smoothed)
    u = u / np.sum(1.0 * u)

    prior_factor = (1 - gamma) * prior_prob_smoothed + gamma * u
    prior_factor = np.power(prior_factor, -alpha)

    # renormalize
    prior_factor = prior_factor / (np.sum(prior_factor * prior_prob_smoothed))

    file_name = os.path.join(data_dir, "CelebA_%s_prior_factor.npy" % size)
    np.save(file_name, prior_factor)

    if do_plot:
        plt.plot(prior_factor)
        plt.yscale("log")
        plt.show()

def bargraph_from_dict(h_dict, group_labels, log=False, colors=None, save=None,
                       title=None, xlabel=None, ylabel=None, legend=False, loc='best'):
    bars = len(h_dict) + 2 # Add one for a space between groups
    groups = len(group_labels)

    print h_dict
    for v in h_dict.values():
        if len(v) != groups:
            raise ValueError("h_dict must contain sequences with the same number as labels")

    plt.cla()
    plts = []
    if log:
        plt.yscale("log")

    if colors == "greyscale":
        colors = [str(i*(.75/len(h_dict))) for i in xrange(len(h_dict))]
    elif colors is None or colors == "spectrum":
        colors = "rygbcmk"

    keys = sorted(h_dict.keys())
    for n, k in enumerate(keys):
        plts.append(plt.bar(range(n+1, bars*groups+1, bars), h_dict[k],
                            color=colors[n % len(colors)]))
    plt.xticks([1 + i*bars + bars/2 for i in xrange(groups)], group_labels)

    if log:
        ymin, ymax = plt.ylim()
        ymin = min(reduce(lambda acc, x: acc+x, h_dict.values(), []))/2.0
        plt.ylim((ymin, ymax))
    if legend:
        plt.legend( [p[0] for p in plts], [k.replace('_', '\_') for k in keys], loc=loc)
    if title:
        plt.title(title)
    if ylabel:
        plt.ylabel(ylabel)
    if xlabel:
        plt.xlabel(xlabel)
    if save is not None:
        plt.savefig(save)

def analyze_neuron(cell, mapping, Mea, signal_range):
    start_t, stop_t = signal_range
    start_t_ixd = np.argmin(np.abs(cell.tvec - start_t))
    stop_t_ixd  = np.argmin(np.abs(cell.tvec - stop_t))
    t_array = cell.tvec[start_t_ixd:stop_t_ixd]
    imem = cell.imem[:, start_t_ixd:stop_t_ixd]
    amps_at_comps = find_amplitude_at_comp(imem, debug=False)
    comp_dist_from_soma = find_comps_dist_from_soma(cell)
    #pl.figure()
    #pl.plot(comp_dist_from_soma, amps_at_comps, 'o')
    #pl.show()


    for elec in xrange(Mea.n_elecs):
        #if not elec == 110:
        #    continue
        print elec
        comp_dist_from_elec = find_comp_dist_from_elec(elec, cell, Mea)
        comp_impact_on_elec = mapping[elec,:]*amps_at_comps
        if 0:
            print np.argmax(comp_dist_from_elec)
            print np.argmin(comp_dist_from_elec)
            pl.figure()
            pl.axis('equal')
            pl.plot(cell.zmid/1000, cell.ymid/1000, 'o')
            pl.plot(Mea.elec_z, Mea.elec_y, 'o')
            pl.plot(cell.zmid[381]/1000, cell.ymid[381]/1000, 'D')
            pl.plot(cell.zmid[791]/1000, cell.ymid[791]/1000, 'D')
            pl.plot(Mea.elec_z[elec], Mea.elec_y[elec], 'x')
            pl.show() 
        pl.plot(comp_dist_from_elec, comp_impact_on_elec, 'o')
        
        pl.yscale('log')
        pl.xscale('log')
        pl.title('Electrode at distal tuft dendrite, when spike originates in tuft,\n'\
             +'and does not propagate to soma.')
        pl.xlabel('Compartment distance from electrode [mm]')
        pl.ylabel('Compartment impact on electrode [uV]')
    pl.show()

def plot_results(list_log, to_plot="losses"):

    list_color = [u'#E24A33',
                  u'#348ABD',
                  u'#FBC15E',
                  u'#777777',
                  u'#988ED5',
                  u'#8EBA42',
                  u'#FFB5B8']

    plt.figure()
    for idx, log in enumerate(list_log):
        with open(log, "r") as f:
            d = json.load(f)

            experiment_name = d["experiment_name"]
            color = list_color[idx]

            plt.plot(d["train_%s" % to_plot],
                     color=color,
                     linewidth=3,
                     label="Train %s" % experiment_name)
            plt.plot(d["val_%s" % to_plot],
                     color=color,
                     linestyle="--",
                     linewidth=3,)
    plt.ylabel(to_plot, fontsize=20)
    if to_plot == "losses":
        plt.yscale("log")
    if to_plot == "accs":
        plt.ylim([0, 1.1])
    plt.xlabel("Number of epochs", fontsize=20)
    plt.title("%s experiment" % dataset, fontsize=22)
    plt.legend(loc="best")
    plt.tight_layout()
    plt.savefig("./figures/%s_results_%s.png" % (dataset, to_plot))
    plt.show()

def plotter(infiles, numcases, labels, linestyles, colors, savefile, ylog=False):
    """Assuming we have output from CORBITS examples,
    creates plots of Transit Probability vs Mean Mutual Inclination.
    """
    plt.figure()
    plt.xlabel('Mean Mutual Inclination [Degrees]', fontsize=16)
    plt.ylabel('Transit Probability', fontsize=16)
    if numcases == 1:
        data = np.loadtxt(infiles[0])
        mmi = data[:,0]
        tp = data[:,1]
        plt.plot(mmi, tp, linestyle=linestyles[0], color=colors[0], linewidth=4)
    else:
        for i in range(numcases):
            data = np.loadtxt(infiles[i])
            mmi = data[:,0]
            tp = data[:,1]
            plt.plot(mmi, tp, linestyle=linestyles[i], color=colors[i], linewidth=4, label=labels[i])
    if numcases != 1:
        plt.legend(loc='best')
    if ylog == True:
        plt.yscale('log')
    plt.savefig(savefile, format='pdf')
    return

def plot_tec(start=None, finish=None, ax=None, x_axis='EPHEMERIS_TIME',
                fmt='k.', bad_fmt='r.', mew=0., labels=True):
    """Quickie that plots TEC in TECU"""

    if ax is None:
        ax = plt.gca()
    else:
        plt.sca(ax)

    if x_axis == 'EPHEMERIS_TIME':
        if start is None:
            start, finish = plt.xlim()
        else:
            plt.xlim(start, finish)

    t = read_tec(start, finish)

    good = t['FLAG'] == True
    plt.plot(t[x_axis][good], t['TEC'][good] / TECU_electrons_per_m2, fmt, mew=mew)
    plt.plot(t[x_axis][~good], t['TEC'][~good] / TECU_electrons_per_m2, bad_fmt, mew=mew)

    plt.xlabel(x_axis)
    plt.ylabel('TEC / TECU')
    plt.yscale('log')