def test_recover_gs_vary_n(Q,baselines,lms):
    """
    This function runs many tests of the linear regression, varying the 
    amount of noise introduced into the data. I hard-coded the global 
    signal strength to be 1 so that it is easy to compare the magnitude
    of the recovered and true global signals and the amplitude of the noise.
    """
    for jj in n.arange(100):
        gs_diff = n.zeros(20,dtype=complex)
        errors = n.zeros(20)
        n_sigs = n.logspace(-3,1,num=20)
        print n_sigs
        for ii,n_sig in enumerate(n_sigs):
            print ii
            gs_true,gs_recov,err = test_recover_gs(Q,baselines,lms,n_sig=n_sig)
            print gs_true
            print gs_recov
            print gs_true.shape
            gs_diff[ii] = gs_recov[0] - gs_true[0] 
            errors[ii] = err
        p.scatter(n_sigs,n.absolute(gs_diff))
        p.scatter(n_sigs,errors,color="red")
        p.xscale('log')
        p.yscale('log')
        p.xlim(1e-4,1e2)
    p.xlabel('Amplitude of noise relative to global signal\n(I.e. true global signal amplitude is 1)')
    #p.ylabel('Recovered global signal (true gs = 1)')
    p.ylabel('Difference between true and recovered global signal')
    #p.show()
    p.savefig('./figures/circle10_Q_pinv_gs_diff_vs_n.pdf')
    p.clf()

def show_table(table_name,ls="none", fmt="o", legend=False, name="m", do_half=0):
	bt = fi.FITS(table_name)[1].read()
	rgpp = (np.unique(bt["rgp_lower"])+np.unique(bt["rgp_upper"]))/2
	nbins = rgpp.size

	plt.xscale("log")
	colours=["purple", "forestgreen", "steelblue", "pink", "darkred", "midnightblue", "gray", "sienna", "olive", "darkviolet"]
	pts = ["o", "D", "x", "^", ">", "<", "1", "s", "*", "+", "."]
	for i,r in enumerate(rgpp):
		sel = (bt["i"]==i)
		snr = 10** ((np.log10(bt["snr_lower"][sel]) + np.log10(bt["snr_upper"][sel]))/2)

		if do_half==1 and i>nbins/2:
			continue
		elif do_half==2 and i<nbins/2:
			continue
		if legend:
			plt.errorbar(snr, bt["%s"%name][i*snr.size:(i*snr.size)+snr.size], bt["err_%s"%name][i*snr.size:(i*snr.size)+snr.size], color=colours[i], ls=ls, fmt=pts[i], lw=2.5, label="$R_{gpp}/R_p = %1.2f-%1.2f$"%(np.unique(bt["rgp_lower"])[i],np.unique(bt["rgp_upper"])[i]))
		else:
			plt.errorbar(snr, bt["%s"%name][i*snr.size:(i*snr.size)+snr.size], bt["err_%s"%name][i*snr.size:(i*snr.size)+snr.size], color=colours[i], ls=ls, fmt=pts[i], lw=2.5)

	plt.xlim(10,300)
	plt.axhline(0, lw=2, color="k")
	
	plt.xlabel("Signal-to-Noise $SNR_w$")
	if name=="m":
		plt.ylim(-0.85,0.05)
		plt.ylabel("Multiplicative Bias $m \equiv (m_1 + m_2)/2$")
	elif name=="alpha":
		plt.ylabel(r"PSF Leakage $\alpha \equiv (\alpha _1 + \alpha _2)/2$")
		plt.ylim(-0.5,2)



	plt.legend(loc="lower right")

def _show_rates(rate, wo, wt, attenuator, tau_NP, tau_P):
    import pylab

    #pylab.figure()
    pylab.errorbar(rate, wt[0], yerr=wt[1], fmt='g.', label='attenuated')
    pylab.errorbar(rate, wo[0], yerr=wo[1], fmt='b.', label='unattenuated')

    pylab.xscale('log')
    pylab.yscale('log')
    pylab.xlabel('incident rate (counts/second)')
    pylab.ylabel('observed rate (counts/second)')
    pylab.legend(loc='best')
    pylab.grid(True)
    pylab.plot(rate, rate/attenuator, 'g-', label='target')
    pylab.plot(rate, rate, 'b-', label='target')

    Ipeak, Rpeak = peak_rate(tau_NP=tau_NP, tau_P=tau_P)
    if rate[0] <= Ipeak <= rate[-1]:
        pylab.axvline(x=Ipeak, ls='--', c='b')
        pylab.text(x=Ipeak, y=0.05, s=' %g'%Ipeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_xaxis_transform())
    if False:
        pylab.axhline(y=Rpeak, ls='--', c='b')
        pylab.text(y=Rpeak, x=0.05, s=' %g\n'%Rpeak,
                   ha='left', va='bottom',
                   transform=pylab.gca().get_yaxis_transform())

 def _set_axis_parameter( self ):
     # set axis to equal length
     params = self.__params
     ax_0 = self._get_axis()
     # set axis aspect 
     pylab.xlim(params['xlim'])
     pylab.ylim(params['ylim'])
     x0,x1 = ax_0.get_xlim()
     y0,y1 = ax_0.get_ylim()
     if params['xlog'] and params['ylog']:
         delta_x = float(np.log(x1)-np.log(x0))
         delta_y = float(np.log(y1)-np.log(y0))
     else:
         delta_x = float(x1 - x0)
         delta_y = float(y1 - y0)
     ax_0.set_aspect(delta_x/delta_y)
     # set tick size
     ax_0.tick_params(axis='both', labelsize=params['ticksize'])
     # set logarithmic scale
     if params['xlog']:
         pylab.xscale('log')
         if params['grid']:
             ax_0.xaxis.grid( True, which='both' )
     if params['ylog']:
         pylab.yscale('log')
         if params['grid']:
             ax_0.yaxis.grid( True, which='both' )
     # grid below bars and boxes
     ax_0.set_axisbelow(params['axisbelow'])

def Validation():
  numSamples = 1000000
  
  theta = np.random.rand(numSamples)*np.pi
  ECo60 = np.array([1.117,1.332])
  Ef0,Ee0 = Compton(ECo60[0],theta)
  Ef1,Ee1 = Compton(ECo60[1],theta)
  dSdE0 = diffXSElectrons(ECo60[0],theta)
  dSdE1 = diffXSElectrons(ECo60[1],theta)

  # Sampling Values
  values = list()
  piMax = np.max([dSdE0,dSdE1])
  while (len(values) < numSamples):
    values.append(SampleRejection(piMax,ComptonScattering))
  # Binning the data
  bins = np.logspace(-3,0.2,100)
  counts = np.histogram(values,bins)
  counts = counts[0]/float(len(values))
  binCenters = 0.5*(bins[1:]+bins[:-1])
  
  # Plotting
  pylab.figure()
  pylab.plot(binCenters,counts,ls='steps')
  #pylab.bar(binCenters,counts,align='center')
  pylab.grid(True)
  pylab.xlim((1E-3,1.4))
  pylab.xlabel('Electron Energy (MeV)')
  pylab.ylabel('Frequency per Photon')
  pylab.yscale('log')
  pylab.xscale('log')
  pylab.savefig('ValComptonScatteringXS.png')

def show_plot(xlabel, ylabel, xlog=False, ylog=False):
  plt.xscale('log') if xlog else None
  plt.yscale('log') if ylog else None
  plt.xlabel(xlabel)
  plt.ylabel(ylabel)
  plt.subplot(111).legend()
  plt.show()

def run_analysis(filename,mode,method):
    click.echo('Reading file : %s'%filename)
    data = IOfile.parsing_input_file(filename)
    click.echo('Creating class...')
    theclass = TFC(data)
    click.echo('Calculating transfer function using %s method'%method)
    if method=='tf_kramer286_sh':
        theclass.tf_kramer286_sh()
    elif method=='tf_knopoff_sh':
        theclass.tf_knopoff_sh()
    elif method=='tf_knopoff_sh_adv':
        theclass.tf_knopoff_sh_adv()
        
    plt.plot(theclass.freq,np.abs(theclass.tf[0]),label=method)
    plt.xlabel('frequency (Hz)')
    plt.ylabel('Amplification')
    plt.yscale('log')
    plt.xscale('log')
    plt.grid(True,which='both')
    plt.legend(loc='best',fancybox=True,framealpha=0.5)
    #plt.axis('tight')
    plt.autoscale(True,axis='x',tight=True)
    plt.tight_layout()
    plt.savefig('test.png', format='png')
    click.echo(click.style('Calculation has been finished!',fg='green'))

def modality(handle, cutoff, resolution):
    """
    """
    x, y = zip(*map(lambda x: map(int, x.split()),
        handle.readlines()[:cutoff]))

    # Ad-hoc weighing function.
    w = map(lambda x: 1.0 / (10 * x + 1), x)
    # Get an approximation of the distribution.
    f = interpolate.UnivariateSpline(x, y, w=w)

    xs = pylab.linspace(0, cutoff, resolution)

    # Plot the original.
    pylab.plot(x, y)

    # Plot the interpolated function.
    ys = f(xs)
    pylab.plot(xs, ys)

    # Plot the tops.
    ys = f(xs, 1)
    g = interpolate.InterpolatedUnivariateSpline(xs, ys)
    pylab.plot(g.roots(), f(g.roots()), 'o')

    # Plot the bending points.
    ys = f(xs, 2)
    g = interpolate.InterpolatedUnivariateSpline(xs, ys)
    pylab.plot(g.roots(), f(g.roots()), 'o')

    pylab.legend(('original', 'interpolated', '1st derivative',
        '2nd derivative'), loc='best')
    pylab.xscale("log")
    pylab.yscale("log")
    pylab.show()

 def plot(self, response, label=None, clip=-80, nharmonics=8, spectrum=None, freq_range=None):
     if freq_range is not None:
         lower_freq, upper_freq = freq_range
     else:
         lower_freq = upper_freq = None
     if lower_freq is None:
         lower_freq = self.start_freq
     if upper_freq is None:
         upper_freq = self.stop_freq
     if self.has_harmonics() and nharmonics > 1 and (spectrum is None or spectrum):
         lines = self.plot_harmonic_spectrum(
             response, nharmonics=nharmonics, lower_freq=lower_freq, upper_freq=upper_freq
         )
         pylab.xlim(left=lower_freq, right=upper_freq)
         for i, line in enumerate(lines):
             if label:
                 line.set_label("%d, %s" % (i + 1, label))
             else:
                 line.set_label("%d" % (i + 1))
         return lines
     elif (spectrum is None and self.has_spectrum()) or spectrum:
         f = numpy.logspace(numpy.log10(lower_freq), numpy.log10(upper_freq), 200)
         h = self.get_spectrum(response, 2 * numpy.pi * f / self.fs)
         s = 20 * numpy.log10(abs(h))
         lines = pylab.plot(f, numpy.where(s > clip, s, numpy.nan), label=label)
         pylab.xscale("log")
         return lines
     else:
         return pylab.plot(self.timeline, response, label=label)

def main():
    # gutenberg
    gu_words = gutenberg.words()
    gu_words_exclude_stops = exclude_stopwords(gu_words)
    gu_fd1 = get_frequency_distribution(gu_words)
    gu_fd2 = get_frequency_distribution(gu_words_exclude_stops)

    pylab.plot(gu_fd1, color='red')
    pylab.plot(gu_fd2, color='orange')

    # inaugural
    in_words = inaugural.words()
    in_words_exclude_stops = exclude_stopwords(in_words)
    in_fd1 = get_frequency_distribution(in_words)
    in_fd2 = get_frequency_distribution(in_words_exclude_stops)

    pylab.plot(in_fd1, color='black')
    pylab.plot(in_fd2, color='gray')

    # reuters
    yen_words = reuters.words(categories='yen')
    yen_words_exclude_stops = exclude_stopwords(yen_words)
    yen_fd1 = get_frequency_distribution(yen_words)
    yen_fd2 = get_frequency_distribution(yen_words_exclude_stops)

    pylab.plot(yen_fd1, color='blue')
    pylab.plot(yen_fd2, color='green')

    pylab.xscale('log')
    pylab.yscale('log')
    pylab.show()

def plot_bins(Y, loglog, logbins, suffix):
    X1, Y1, X2, Y2, alpha = fit_data(Y, True, pdf=True)      

    
    
    pylab.figure(figsize=(7.5, 7))
    
    pylab.rcParams.update({'font.size': 20})

    pylab.scatter(X1, Y1)

    pylab.plot(X2, Y2, '--')

    bounds = get_bounds(X1, Y1, loglog, loglog)

    if loglog:
        pylab.xscale('log')
        pylab.yscale('log')
        xtext = numpy.exp(numpy.log(bounds[0])+(numpy.log(bounds[1])-numpy.log(bounds[0]))*0.65)
        ytext = numpy.exp(numpy.log(bounds[2])+(numpy.log(bounds[3])-numpy.log(bounds[2]))*0.65)
    else:
        xtext = (bounds[0]+bounds[1])/2.0
        ytext = (bounds[2]+bounds[3])/2.0

    pylab.axis(bounds)

    pylab.text(xtext, ytext, '$gamma$='+'{0:.2f}'.format(alpha))
    
    pylab.xlabel('Change')
    pylab.ylabel('Density')

    pylab.tight_layout()

    pylab.show()

    def plot(self, ylog10scale = False, timescale = "years", year = 25):
        """
        Generate figure and axis for the population structure
        timescale choose from "2N0", "4N0", "generation" or "years"
        """
        time = self.Time
        pop  = self.pop
        for i in range(1,len(self.pop)):
            if type(pop[i]) == type(""):
                # ignore migration commands, and replace by (unchanged) pop size
                pop[i] = pop[i-1]
        
        if time[0] != 0 :
            time.insert(0, float(0))
            pop.insert(0, float(1))
        
        if timescale == "years":
            time = [ti * 4 * self.scaling_N0 * year for ti in time ]
            pl.xlabel("Time (years, "+`year`+" years per generation)",  fontsize=20)    
            #pl.xlabel("Years")    
        elif timescale == "generation":
            time = [ti * 4 * self.scaling_N0 for ti in time ]
            pl.xlabel("Generations)")    
        elif timescale == "4N0":
            time = [ti*1 for ti in time ]
            pl.xlabel("Time (4N generations)")    
        elif timescale == "2N0":
            time = [ti*2 for ti in time ]
            pl.xlabel("Time (2N generations)")       
        else:
            print "timescale must be one of \"4N0\", \"generation\", or \"years\""
            return
        
        time[0] = time[1] / float(20)
        
        time.append(time[-1] * 2)
        yaxis_scaler = 10000
        
        pop = [popi * self.scaling_N0 / float(yaxis_scaler) for popi in pop ]
        pop.insert(0, pop[0])               
        pl.xscale ('log', basex = 10)        
        #pl.xlim(min(time), max(time))
        pl.xlim(1e3, 1e7)
        
        if ylog10scale:
            pl.ylim(0.06, 10000)
            pl.yscale ('log', basey = 10)            
        else:
            pl.ylim(0, max(pop)+2)
        
        pl.ylim(0,5)            
        pl.tick_params(labelsize=20)

        #pl.step(time, pop , color = "blue", linewidth=5.0)
        pl.step(time, pop , color = "red", linewidth=5.0)
        pl.grid()
        #pl.step(time, pop , color = "black", linewidth=5.0)
        #pl.title ( self.case + " population structure" )
        #pl.ylabel("Pop size ($*$ "+`yaxis_scaler` +")")
        pl.ylabel("Effective population size",fontsize=20 )

def plot_fas(freqs, ns_data, ew_data, eas_smoothed_data, fas_plot, station):
    """
    Create a plot of both FAS components
    """
    # Generate plot

    # Set plot dims
    pylab.gcf().set_size_inches(11, 8.5)
    pylab.gcf().clf()

    # Adjust title y-position
    t = pylab.title("Station: %s" % (station), size=12)

    pylab.plot(freqs, ns_data, 'b', lw=0.75, label="NS")
    pylab.plot(freqs, ew_data, 'r', lw=0.75, label="EW")
    pylab.plot(freqs, eas_smoothed_data, 'k', lw=1.25, label="Smoothed EAS")
    pylab.legend(loc='upper right')
    pylab.xscale('log')
    pylab.yscale('log')
    pylab.ylabel('Fourier Amplitude (cm/s)')
    pylab.xlabel('Frequency (Hz)')
    pylab.axis([0.01, 100, 0.001, 1000])
    pylab.grid(True)
    pylab.grid(b=True, which='major', linestyle='-', color='lightgray')
    pylab.grid(b=True, which='minor', linewidth=0.5, color='gray')

    # Save plot
    pylab.savefig(fas_plot, format="png",
                  transparent=False, dpi=plot_config.dpi)
    pylab.close()

    def test_fit_bb(self):
        def func(nu, T):
            return np.pi * rf.planck(nu, T, inp="Hz", out="freq")

        self.sp.cut_flux(max(self.sp.Flux) * 1e-5)

        freq = self.sp.Freq
        flux = self.sp.Flux
        Tinit = 1.e4

        popt, pcov = curve_fit(func, freq, flux, p0=Tinit)
        Tbest = popt
        # bestT, pcov = curve_fit(rf.fit_planck(nu, inp='Hz'), nu, flux, p0=Tinit, sigma=sigma)
        sigmaT = np.sqrt(np.diag(pcov))

        print 'True model values'
        print '  Tbb = %.2f K' % self.sp.T

        print 'Parameters of best-fitting model:'
        print '  T = %.2f +/- %.2f K' % (Tbest, sigmaT)

        # Tcol = self.sp.temp_color
        ybest = np.pi * rf.planck(freq, Tbest, inp="Hz", out="freq")
        # plot the solution
        plt.plot(freq, flux, 'b*', label='Spectral T: %f' % self.sp.T)
        plt.plot(freq, ybest, 'r-', label='Best Tcol: %f' % Tbest)
        plt.xscale('log')
        plt.yscale('log')
        plt.legend(loc=3)
        plt.show()

        self.assertAlmostEqual(Tbest, self.sp.T,
                               msg="For planck  Tcolor [%f] should be equal sp.T [%f]." % (Tbest, self.sp.T),
                               delta=Tbest*0.01)

def filterresponse(b,a,fs=44100,scale='log',**kwargs):
    w, h = freqz(b,a)
    pl.subplot(2,1,1)
    pl.title('Digital filter frequency response')
    pl.plot(w/max(w)*fs/2, 20 * np.log10(abs(h)),**kwargs)
    pl.xscale(scale)
#    if scale=='log':
#        pl.semilogx(w/max(w)*fs/2, 20*np.log10(np.abs(h)), 'k')
#    else:
#        pl.plot(w/max(w)*fs/2, 20*np.log10(np.abs(h)), 'k')
        
    pl.ylabel('Gain (dB)')
    pl.xlabel('Frequency (rad/sample)')
    pl.axis('tight')    
    pl.grid()    
    

    pl.subplot(2,1,2)
    angles = np.unwrap(np.angle(h))
    if scale=='log':
        pl.semilogx(w/max(w)*fs/2, angles, **kwargs)
    else:
        pl.plot(w/max(w)*fs/2, angles, **kwargs)
        

    pl.ylabel('Angle (radians)')
    pl.grid()
    pl.axis('tight')
    pl.xlabel('Frequency (rad/sample)')

def plot_values(X, Y, xlabel, ylabel, suffix):
    output_filename = constants.CHARTS_FOLDER_NAME + constants.DATASET + '_' + suffix

    
    pylab.figure(figsize=(8, 7))

    pylab.rcParams.update({'font.size': 20})

    pylab.scatter(X, Y)
    
    '''
    #smoothing
    s = np.square(np.max(Y))
    tck = interpolate.splrep(X, Y, s=s)
    Y_smooth = interpolate.splev(X, tck)
    pylab.plot(X, Y_smooth)
    '''

    #pylab.axis(vis.get_bounds(X, Y, False, False))

    pylab.xscale('log')
    pylab.yscale('log')

    pylab.xlabel(xlabel)
    pylab.ylabel(ylabel)   

    #pylab.tight_layout()

    pylab.savefig(output_filename + '.pdf')

    def mesh2d_mcolor_mask(self, data, axis, output=None, mask=None, datscale='log', 
           axiscale=['log', 'log'], pcolors='Greys', maskcolors=None):

        """       >>> generate 2D mesh plot <<<
	"""

        pl.clf()
	fig=pl.figure()
	ax=fig.add_subplot(111)

	pldat=data  
	
        # get the color norm
	if(datscale=='log'):
	    cnorm=colors.LogNorm()
	elif(datscale=='linear'):
	    cnorm=colors.NoNorm()
	else:
	    raise Exception


        color1=colors.colorConverter.to_rgba('white')
        color2=colors.colorConverter.to_rgba('blue')
        color3=colors.colorConverter.to_rgba('yellow')
        my_cmap0=colors.LinearSegmentedColormap.from_list('mycmap0',[color1, color1, color2, color2, color2, color3, color3], 512) 
        my_cmap0._init()


        if pcolors!=None:
            cm=ax.pcolormesh(axis[0,:], axis[1,:], pldat, cmap=pl.cm.get_cmap(pcolors),
	                     norm=cnorm) 

            #cm=ax.pcolormesh(axis[0,:], axis[1,:], pldat, cmap=my_cmap0, norm=cnorm) 
	else:
            cm=ax.pcolormesh(axis[0,:], axis[1,:], pldat, norm=cnorm) 


        if mask!=None:

            # get the color map of mask
	    """
            color1=colors.colorConverter.to_rgba('white')
            color2=colors.colorConverter.to_rgba('red')
            my_cmap=colors.LinearSegmentedColormap.from_list('mycmap',[color1, color2], 512) 
            my_cmap._init()
            alphas=np.linspace(0.2, 0.7, my_cmap.N+3)
            my_cmap._lut[:,-1] = alphas 
	    """
    
	    maskdata=np.ma.masked_where((mask<=1e-2)&(mask>=-1e-2) , mask)
            mymap=ax.contourf(axis[0,:], axis[1,:], maskdata, cmap=maskcolors)

            cbar=fig.colorbar(mymap, ticks=[4, 6, 8]) #, orientation='horizontal')
            cbar.ax.set_yticklabels(['void', 'filament', 'halo'])

	pl.xscale(axiscale[0])
	pl.yscale(axiscale[1])


        return 

def plot_xy(cursor, query, prefix=None, color='b', marker='.', xlog=False, ylog=False, xlabel='', ylabel='', title=''):
    """
        Executes the 'query' which should return two numerical columns.
    """
    cursor.execute(query)
    x_list = []
    y_list = []
    for row in cursor:
        (x, y) = row
        if (x != None and y != None):
            x_list.append(x)
            y_list.append(y)
    
    X = pylab.array(x_list)
    Y = pylab.array(y_list)
    pylab.figure()
    pylab.hold(True)
    pylab.plot(X, Y, color=color, marker=marker, linestyle='None')
    if (xlog):
        pylab.xscale('log')
    if (ylog):
        pylab.yscale('log')
    
    pylab.title(title + " (R^2 = %.2f)" % pylab.corrcoef(X,Y)[0,1]**2)
    pylab.xlabel(xlabel)
    pylab.ylabel(ylabel)
    if (prefix != None):
        pylab.savefig('../res/%s.pdf' % prefix, format='pdf')
    pylab.hold(False)

def test_plot_schechter():

    phiStar = 1.8e-3
    MStar = -20.04
    alpha = -1.71
    
    LStar = magnitudes.L_nu_from_magAB(MStar)

    mags = numpy.arange(-22, -11.0, 0.5)
    lums = magnitudes.L_nu_from_magAB(mags)

    phi_L = schechterL(lums, phiStar, alpha, LStar)
    phi_M = schechterM(mags, phiStar, alpha, MStar)

    L_L = schechterCumuLL(lums, phiStar, alpha, LStar)
    L_M = schechterCumuLM(mags, phiStar, alpha, MStar)

    phi_L_func = lambda l: l * schechterL(l, phiStar, alpha, LStar)
    L_L_num = utils.logquad(phi_L_func, lums, 1e35)[0]
    L_L_num2 = utils.vecquad(phi_L_func, lums, 1e29)[0]
    
    phi_M_func = lambda m: (magnitudes.L_nu_from_magAB(m) *
                            schechterM(m, phiStar, alpha, MStar))
    L_M_num2 = utils.vecquad(phi_M_func, -25, mags)[0]

    Ltot_L = schechterTotLL(phiStar, alpha, LStar)
    Ltot_M = schechterTotLM(phiStar, alpha, MStar)

    pylab.figure()
    pylab.subplot(221)
    pylab.plot(lums, lums * lums * phi_L)
    pylab.xscale('log')
    pylab.yscale('log')
    pylab.ylabel(r'$ L^2 \Phi_L$')

    pylab.subplot(222)
    pylab.plot(mags, -mags * lums * phi_M)
    pylab.yscale('log')
    pylab.ylabel(r'$ -M L \Phi_M$')

    pylab.subplot(223)
    pylab.plot(lums, Ltot_L - L_L)
    pylab.plot(lums, L_M)
    pylab.plot(lums, L_L_num, '--')
    pylab.plot(lums, L_L_num2, ':')
    pylab.plot(lums, L_M_num2, 'x')
    pylab.axhline(y=Ltot_L)
    pylab.axhline(y=Ltot_M)
    pylab.xscale('log')
    pylab.yscale('log')

    pylab.subplot(224)
    pylab.plot(mags, Ltot_M - L_M)
    pylab.plot(mags, L_L)
    pylab.plot(mags, L_L_num, '--')
    pylab.plot(mags, L_L_num2, ':')
    pylab.plot(mags, L_M_num2, 'x')
    pylab.axhline(y=Ltot_L)
    pylab.axhline(y=Ltot_M)
    pylab.yscale('log')

def plot_per_sel(performance, selected, top, repeats, xlabel, ylabel, suffix):
    print suffix

    output_filename = constants.WISDOM_FOLDER_NAME + suffix

    pylab.figure(figsize=(15, 10))
    
    pylab.rcParams.update({'font.size': 30})

    X, Y = aggregate_per_sel(performance, selected, top, repeats, 'true', '0')
    pylab.plot(X, Y, linewidth=2, color="red", marker='o', markersize=15, markeredgecolor="red", markerfacecolor="white")

    X, Y = aggregate_per_sel(performance, selected, top, repeats, 'true', '1')
    pylab.plot(X, Y, linewidth=2, color="blue", marker='o', markersize=15, markeredgecolor="blue", markerfacecolor="white")

    X, Y = aggregate_per_sel(performance, selected, top, repeats, 'true', '-1')
    pylab.plot(X, Y, linewidth=2, color="green", marker='o', markersize=15, markeredgecolor="green", markerfacecolor="white")

    pylab.xscale('log', basex=10)

    pylab.xlabel(xlabel)
    pylab.ylabel(ylabel)

    pylab.legend(['Crowd', 'Network - 1', 'Network - Inf'], loc='upper left')
    
    pylab.tight_layout()

    pylab.savefig(output_filename + '.pdf')