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_degreeRate(db, keynames, save_path):
	degRate_x_name = 'degRateDistr_x'
	degRate_y_name = 'degRateDistr_y'

	plt.clf()
	plt.figure(figsize = (8, 5))

	plt.subplot(1, 2, 1)
	plt.plot(db[keynames['mog']][degRate_x_name], db[keynames['mog']][degRate_y_name], 'b-', lw=5, label = 'fairyland')
	plt.plot(db[keynames['mblg']][degRate_x_name], db[keynames['mblg']][degRate_y_name], 'r:', lw=5, label = 'twitter')
	plt.plot(db[keynames['im']][degRate_x_name], db[keynames['im']][degRate_y_name], 'k--', lw=5, label = 'yahoo')
	plt.xscale('log')
	plt.grid(True)
	plt.title('interaction')
	plt.legend(('fairyland', 'twitter', 'yahoo'), loc = 4, prop = {'size': 10})
	plt.xlabel('In-degree to Out-degree Ratio')
	plt.ylabel('CDF')

	plt.subplot(1, 2, 2)
	plt.plot(db[keynames['mogF']][degRate_x_name], db[keynames['mogF']][degRate_y_name], 'b-', lw=5, label = 'fairyland')
	plt.plot(db[keynames['mblgF']][degRate_x_name], db[keynames['mblgF']][degRate_y_name], 'r:', lw=5, label = 'twitter')
	#plt.plot(db[keynames['imF']][degRate_x_name], db[keynames['imF']][degRate_y_name], 'k--', lw=5, label = 'yahoo')
	plt.xscale('log')
	plt.grid(True)
	plt.title('ally')
	plt.xlabel('In-degree to Out-degree Ratio')
	plt.ylabel('CDF')

	plt.savefig(os.path.join(save_dir, save_path))

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 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 wykres(katalog):
    import numpy as np
    import os
    import scipy.io.wavfile as siw
    import json
    import scipy.fftpack as sf
    import matplotlib.pylab as plt
    os.chdir(katalog)
    song = open('song.txt', 'r')
    songlist = list()
    defs = json.load(open('defs.txt', 'r'))
    beat = defs["bpm"]
    nframes = 60/beat*44100
    for songline in song:
        track = open('track'+songline.strip('\n')+'.txt', 'r')
        tracklist = list()
        for trackline in track:
            y = np.zeros([nframes, 2])
            for i in range(len(trackline.split())):
                y = y + siw.read('sample'+''.join(trackline.split()[i])+'.wav')[1][0:nframes]
            y = np.mean(y,axis=1)
            y /= 32767
            tracklist = np.r_[tracklist, y]
        songlist = np.r_[songlist, tracklist]
    yf = sf.fft(songlist)
    n = yf.size
    xf = np.linspace(0, 44100/2, n/2)
    plt.plot(xf, 2*(yf[1:int(n/2)+1])/n)
    plt.xscale('log')
    plt.show()

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 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 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 plot_scatter(times):
    pl.scatter(times[:, 0], times[:, 1])
    pl.xlabel('Time in seconds')
    pl.ylabel('Ratio')
    pl.title('Time in seconds versus the ratio')
    pl.xlim(0, 300)
    pl.xscale('log')
    pl.show()

def plot_times():
    res = []
    with open('times', 'r') as f:
        for l in f.readlines():
            x = eval(l)
            res += [(x['samples'], x['auto'], x['tri'])]

    res = list(zip(*res))
    pl_auto = pl.plot(res[0], res[1])
    pl_tri = pl.plot(res[0], res[2])
    pl.xscale('log')
    pl.ylabel("time (in seconds)")
    pl.xlabel("points (logarithmic scale)")
    pl.legend((pl_auto[0], pl_tri[0]), ('auto', 'tri'))

    pl.show()

def doPlot(parallel_aucs, serial_aucs, times, errors):
    times = list(times)
    times_histo = np.histogram(parallel_aucs,bins=times)
    #values,edges = times_histo
    parallel_values = parallel_aucs[1:]
    edges = times
    print(len(parallel_values))
    print(len(edges))
    serial_values = np.array(serial_aucs[1:])
    errors = np.array(errors[1:])
    edges = np.array(times[:-1])
    print(errors.shape)
    print(edges.shape)
    print(serial_values.shape)


    plt.figure()
    plt.plot(edges, parallel_values,label = "Distributed search") #, width=np.diff(edges), ec="k", align="edge")
    plt.plot(edges, serial_values, label="Sequential search") #, width=np.diff(edges), ec="k", align="edge")
    #plt.fill_between(edges, serial_values-errors,serial_values+errors)
    plt.legend(loc = (0.6,0.7))
    plt.xlabel("Time [minutes]", fontsize=20)
    #plt.yscale('log')
    plt.ylabel('Best validation AUC', fontsize=20)
    plt.savefig("times.png")

    plt.figure()
    plt.plot(edges, parallel_values,label = "Distributed search") #, width=np.diff(edges), ec="k", align="edge")
    plt.plot(edges, serial_values, label="Sequential search") #, width=np.diff(edges), ec="k", align="edge")
    #plt.fill_between(edges, serial_values-errors,serial_values+errors)
    plt.legend(loc = (0.6,0.7))
    plt.xlabel("Time [minutes]", fontsize=20)
    plt.xscale('log')
    plt.xlim([0,100])
    plt.ylabel('Best validation AUC', fontsize=20)
    plt.savefig("times_logx_start.png")

    plt.figure()
    plt.plot(edges, parallel_values,label = "Distributed search") #, width=np.diff(edges), ec="k", align="edge")
    plt.plot(edges, serial_values, label="Sequential search") #, width=np.diff(edges), ec="k", align="edge")
    #plt.fill_between(edges, serial_values-errors,serial_values+errors)
    plt.legend(loc = (0.6,0.7))
    plt.xlabel("Time [minutes]", fontsize=20)
    plt.xscale('log')
    plt.xlim([100,10000])
    plt.ylabel('Best validation AUC', fontsize=20)
    plt.savefig("times_logx.png")

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 plot(self,bins=1000,der=0,log=False,error=False,color=None):
        if log==True:
            x=np.logspace(log10(self.minX)+1e-7,log10(self.maxX)-1e-7,num=bins)
            plt.xscale("log")
        else:
            x=np.linspace(self.minX,self.maxX,num=bins)

        y=self.__call__(x,der=der)
        plt.plot(x,y,color=color)
        
        if(error==True):
            ymin=self.__call__(x,der=der,value="top")
            ymax=self.__call__(x,der=der,value="bottom")
            
            plt.plot(x,ymax,color=color)
            plt.plot(x,ymin,color=color)
            plt.fill_between(x, ymin,ymax,alpha=0.5,color=color)

def main():
    arg_parser = argparse.ArgumentParser(
        description="extracts statistical information from a given CSV video file created with video2csv"
    )
    arg_parser.add_argument("--csvFile", help="the movie file to convert with ffmpeg")
    args = arg_parser.parse_args()

    if args.csvFile is None:
        print("no CSV file specified!")
    if not os.path.isfile(args.csvFile):
        print("No CSV file given or file does not exist")
        sys.exit(127)

    iframes = []
    pframes = []

    with open(args.csvFile) as csvfile:
        video_reader = csv.reader(csvfile, delimiter=";")
        video_reader.next()  # skip first line comment
        for row in video_reader:
            if "I" in row[1]:

                iframes.append(float(row[3]))
            else:

                pframes.append(float(row[3]))

        print("Number of I-frames: %s" % (len(iframes)))
        print("Average I-frame size: %s" % (np.average(iframes)))
        print("Number of P-frames: %s" % (len(pframes)))
        print("Average P-frame size: %s" % (np.average(pframes)))

        x1, y1 = list_to_ccdf(iframes)
        x2, y2 = list_to_ccdf(pframes)
        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)

        ax.plot(x1, y1, label="IFrames")
        ax.plot(x2, y2, label="PFrames")
        plt.xscale("log")
        plt.xlabel("coded picture size in bytes")
        plt.ylabel(("P"))
        plt.grid()
        plt.legend()
        plt.show()

def SMBH_mass(save=False):

	fig=plot.figure()

	#a=ascii.read('data/BHmass_dist.dat',names=['mass','N'],data_start=1)
	#mass=(np.round(a['mass']*100.)/100.)
	#N=np.array(np.round(a['N']*100),dtype=np.int64)

	high=ascii.read('data/BH_mass_High_z.txt',names=['mass'])
	low=ascii.read('data/BH_mass_Low_z.txt',names=['mass'])

	loghigh=np.log10(high['mass'])
	loglow=np.log10(low['mass'])
	#ind1=np.arange(0,13,1)
	#ind2=np.arange(13,len(mass),1)

	#m1=np.repeat(mass[ind1],N[ind1])
	#w1=np.repeat(np.repeat(1./max(N[ind1]),len(ind1)),N[ind1])
	#m2=np.repeat(mass[ind2],N[ind2])
	#w2=np.repeat(np.repeat(1./max(N[ind2]),len(ind2)),N[ind2])

	low_bin=np.logspace(np.min(loglow),np.max(loglow),num=14)
	plot.hist(low['mass'],bins=low_bin,color='blue',weights=np.repeat(1./28,len(low)))
	high_bin=np.logspace(np.min(loghigh),np.max(loghigh),num=10)
	plot.hist(high['mass'],bins=high_bin,color='red',alpha=0.7,weights=np.repeat(1./28,len(high)))

	plot.annotate('Low Redshift (z < 3) Blazars',xy=(1.5e9,0.8),xycoords='data',fontsize=14,color='blue')
	plot.annotate('High Redshift (z > 3) Blazars',xy=(1.5e9,0.5),xycoords='data',fontsize=14,color='red')

	plot.xlim([5e7,5e10])
	plot.ylim([0,1.05])
	plot.xscale('log')
#	plot.yscale('log')
	plot.xlabel(r'Black Hole Mass (M$_{\odot}$)')
	plot.ylabel('Fraction of Known Blazars')
#	plot.title('Supermassive Black Hole Mass Evolution')

	if save:
		plot.savefig('SMBH_mass.png', bbox_inches='tight')
		plot.savefig('SMBH_mass.eps', bbox_inches='tight')
	else:
		plot.show()

	return

def plot(result):

    from matplotlib import pylab as pl
    import scipy as sp

    if not result:
        result = sp.loadtxt('morph.csv', delimiter=',', skiprows=1).T

    x, y = result[0], result[1:]

    for i in y:
        pl.plot(x, i)

    pl.xlabel('Number of characters')
    pl.ylabel('Time (sec)')
    pl.xscale('log')
    pl.grid(True)
    pl.savefig("images/time.png")
    pl.show()

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 compute_weights(correlation, verbose=0, uselatex=False):
        """ computes and returns weights of the same length as 
        `correlation.correlation_fit`
        
        `correlation` is an instance of Correlation
        """
        corr = correlation
        model = corr.fit_model
        model_parms = corr.fit_parameters
        ival = corr.fit_ival
        weight_data = corr.fit_weight_data
        weight_type = corr.fit_weight_type
        #parameters = corr.fit_parameters
        #parameters_range = corr.fit_parameters_range
        #parameters_variable = corr.fit_parameters_variable
        
        cdat = corr.correlation
        if cdat is None:
            raise ValueError("Cannot compute weights; No correlation given!")
        cdatfit = corr.correlation_fit
        x_full = cdat[:,0]
        y_full = cdat[:,1]
        x_fit = cdatfit[:,0]
        #y_fit = cdatfit[:,1]
        
        dataweights = np.ones_like(x_fit)

        try:
            weight_spread = int(weight_data)
        except:
            if verbose > 1:
                warnings.warn("Could not get weight spread for spline. Setting it to 3.")
            weight_spread = 3

        if weight_type[:6] == "spline":
            # Number of knots to use for spline
            try:
                knotnumber = int(weight_type[6:])
            except:
                if verbose > 1:
                    print("Could not get knot number. Setting it to 5.")
                knotnumber = 5
            
            # Compute borders for spline fit.
            if ival[0] < weight_spread:
                # optimal case
                pmin = ival[0]
            else:
                # non-optimal case
                # we need to cut pmin
                pmin = weight_spread
            if x_full.shape[0] - ival[1] < weight_spread:
                # optimal case
                pmax = x_full.shape[0] - ival[1]
            else:
                # non-optimal case
                # we need to cut pmax
                pmax = weight_spread

            x = x_full[ival[0]-pmin:ival[1]+pmax]
            y = y_full[ival[0]-pmin:ival[1]+pmax]
            # we are fitting knots on a base 10 logarithmic scale.
            xs = np.log10(x)
            knots = np.linspace(xs[1], xs[-1], knotnumber+2)[1:-1]
            try:
                tck = spintp.splrep(xs, y, s=0, k=3, t=knots, task=-1)
                ys = spintp.splev(xs, tck, der=0)
            except:
                if verbose > 0:
                    raise ValueError("Could not find spline fit with "+\
                                     "{} knots.".format(knotnumber))
                return
            if verbose > 0:
                try:
                    # If plotting module is available:
                    name = "spline fit: "+str(knotnumber)+" knots"
                    plotting.savePlotSingle(name, 1*x, 1*y, 1*ys,
                                             dirname=".",
                                             uselatex=uselatex)
                except:
                    # use matplotlib.pylab
                    try:
                        from matplotlib import pylab as plt
                        plt.xscale("log")
                        plt.plot(x, ys, x, y)
                        plt.show()
                    except ImportError:
                        # Tell the user to install matplotlib
                        print("Couldn't import pylab! - not Plotting")

            ## Calculation of variance
            # In some cases, the actual cropping interval from ival[0]
            # to ival[1] is chosen, such that the dataweights must be
            # calculated from unknown datapoints.
            # (e.g. points+endcrop > len(correlation)
            # We deal with this by multiplying dataweights with a factor
            # corresponding to the missed points.
            for i in range(x_fit.shape[0]):
                # Define start and end positions of the sections from
                # where we wish to calculate the dataweights.
#.........这里部分代码省略.........

     
     # -----------set up units properly------------------#        
     mag_order=np.int((1)*np.round(np.log10(np.mean(sp.data))))
     sp.xarr.units='angstroms'
     sp.xarr.xtype = 'wavelength'
     sp.units = r'$10^{'+str(mag_order)+'}$ erg s$^{-1}$ $cm^{-1}$'
     sp.data /= 10**(mag_order)    
     #-------------- set up units properly------------#
     
     
     wlmin=sp.xarr[0]
     wlmax=sp.xarr[-1]
 
 
     pylab.yscale('log')
     pylab.xscale('log')
     
     #pylab.xscale('log',subsx=[3,4])
     pylab.rcParams["figure.figsize"]=16,6
     sp.plotter(figure=1,xmin=wlmin,xmax=wlmax,ymin=1.1*sp.data.min(),ymax=3.5*sp.data.max())
     pylab.close()
     
     #-------------continuous subtraction-------------------#
     if w=="UV": 
         continuous,wlmin_UV=continuous_substraction( i, sp, mag_order,UV_limits,w)
         backup=sp.copy()
         sp.data=sp.data - continuous
         balmer_template,balmer_tot, index=balmer_normalization(sp,balmer_cont_template,balmer_lines_template)
         
         sp.data=backup.data
         

i = 1
for line in f:
    s = line

    word = s.split()  # separating words and frequencies
    wordList.append(word[0])  # storing words
    freqList.append(word[1])  # storing frequencies
    rankList.append(i)  # storing ranks
    i = i + 1
for i in range(0, len(freqList)):
    freqList[i] = int(freqList[i])
rankList = np.log(rankList)  # converting to log values
freqList = np.log(freqList)

pl.yscale("log")  # defining log-log scale of graph
pl.xscale("log")

pl.plot(rankList, freqList)

# pl.axis([0,75000,164438,23135851162])
po = np.polyfit(rankList, freqList, 5)  # for the smoothed graph
yfit = np.polyval(po, rankList)

pl.plot(rankList, yfit)  # plotting the graph
pl.xlabel("Log Rank")
pl.ylabel("Log Frequency")  # defining labels, titles and legends

pl.suptitle("Solution 3 : Frequency vs Rank")
pl.legend(["Log Frequency vs Log Rank", "Smoothened Graph"])
pl.show()