Python pylab.step函数代码示例

def measure_flexure_y(fine, HDUlist, profwidth=5, plot=False, outname=None):

    dat = HDUlist[0].data
    exptime = HDUlist[0].header["EXPTIME"]

    profs = []
    xx = np.arange(profwidth * 2)

    for ix in np.arange(500, 1200, 10):
        f = fine[ix]
        profile = np.zeros(profwidth * 2)

        # bad fit
        if not f.ok:
            continue
        # noisy fit
        if f.lamnrms > 1:
            continue
        # short spectrum
        if f.xrange[1] - f.xrange[0] < 200:
            continue

        yfun = np.poly1d(f.poly)
        for xpos in np.arange(f.xrange[0], f.xrange[1]):
            try:
                ypos = int(np.round(yfun(xpos)))
            except:
                continue
            try:
                profile += dat[ypos - profwidth : ypos + profwidth, xpos]
            except:
                continue

        profs.append(profile - np.min(profile))

    if plot:
        pl.figure(1)
    ffun = FF.mpfit_residuals(FF.gaussian4)
    parinfo = [{"value": 1}, {"value": profwidth}, {"value": 2}, {"value": 0}, {"value": 0}]

    profposys = []
    for prof in profs:
        if plot:
            pl.step(xx, prof)
        parinfo[0]["value"] = np.max(prof)
        fit = FF.mpfit_do(ffun, xx, prof, parinfo)
        if plot:
            pl.plot(xx, FF.gaussian4(fit.params, xx))
        profposys.append(fit.params[1] - profwidth - 1)
    if plot:
        pl.show()
    profposys = np.array(profposys)

    mn = np.mean(profposys)
    sd = np.std(profposys)
    ok = np.abs(profposys - mn) / sd < 3
    required_shift = np.mean(profposys[ok])
    print "dY = %3.2f pixel shift" % required_shift

    return required_shift

def test_metropolis(Nsamp=int(1e4), Nburnin=0, sampleSig=5, thetaInit=-5):
    from scipy import stats
    import pylab as plt

    def testfunc(theta):
        return (stats.cauchy(-10, 2).pdf(theta)
                + 4 * stats.cauchy(10, 4).pdf(theta))

    def lnp(theta):
        """ interface to metrop """
        return (0, np.log10(testfunc(theta)))

    x = np.linspace(-50, 50, 1e4)
    y = testfunc(x)

    # compute normalization function: used later to put on same as histogram
    Zfunc = np.sum(y * x.ptp() / len(x))

    plt.plot(x, y / Zfunc, 'k-', lw=2)
    plt.xlabel('theta')
    plt.ylabel('f')

    lnps, samples = metropolis(lnp, [thetaInit], Nburnin=Nburnin, Nsamp=Nsamp,
                               sampleCov=sampleSig ** 2)

    n, b = np.histogram(samples.ravel(), bins=np.arange(-50, 50, 1),
                        normed=False)
    Zhist = (n * np.diff(b)).sum()
    plt.step(b[1:], n.astype(float) / Zhist, color='r')
    plt.title('test metropolis')

def plot_attenuation_sim(rmArray, f, days_in_season=82, n_iter=1000, label=''):
	N = days_in_season
	l2 = (0.3 / f) ** 2
	hist, bins = np.histogram(rm, bins=np.arange(min(rmArray), max(rmArray) + bin_width, bin_width))
	bin_midpoints = bins[:-1] + np.diff(bins) / 2
	cdf = np.cumsum(hist)
	cdf = cdf / cdf[-1]
	values = np.random.rand(n_iter)
	value_bins = np.searchsorted(cdf, values)
	random_from_cdf = bin_midpoints[value_bins]

	factors = []

	for i, RM in enumerate(random_from_cdf):
		atten = 0
		for j in range(i + 1, N, 1):
			atten += np.cos(2 * (random_from_cdf[j] - RM) * l2)
		factor = (float(N) + (2. * atten)) / pow(float(N), 2.)
		factors.append(factor)

	Pv, v = np.histogram(factors, bins=10, density=True)
	print('Epsilon (sim) =', np.mean(factors), '+/-', np.std(factors))
	v = 0.5 * (v[1:] + v[:-1])
	pylab.step(v, Pv, label=label)

	return [np.mean(factors), np.std(factors)]

def test_sparse_fused_lasso(n=500,l1=2.,ratio=1.,control=control):

    D = (np.identity(n) - np.diag(np.ones(n-1),-1))[1:]
    D = np.vstack([D, ratio*np.identity(n)])
    M = np.linalg.eigvalsh(np.dot(D.T, D)).max() 

    Y = np.random.standard_normal(n)
    Y[int(0.1*n):int(0.3*n)] += 3.
    p1 = sapprox.signal_approximator((D, Y),L=M)
    p1.assign_penalty(l1=l1)
    
    t1 = time.time()
    opt1 = regreg.FISTA(p1)
    opt1.fit(tol=control['tol'], max_its=control['max_its'])
    beta1 = opt1.problem.coefs
    t2 = time.time()
    ts1 = t2-t1

    beta, _ = opt1.output
    X = np.arange(n)
    if control['plot']:
        pylab.clf()
        pylab.step(X, beta, linewidth=3, c='red')
        pylab.scatter(X, Y)
        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 createPlot(dataY, dataX, ticksX, annotations, axisY, axisX, dostep, doannotate):
    if not ticksX:
        ticksX = dataX
    
    if dostep:
        py.step(dataX, dataY, where='post', linestyle='-', label=axisY) # where=post steps after point
    else:
        py.plot(dataX, dataY, marker='o', ms=5.0, linestyle='-', label=axisY)
    
    if annotations and doannotate:
        for note, x, y in zip(annotations, dataX, dataY):
            py.annotate(note, (x, y), xytext=(2,2), xycoords='data', textcoords='offset points')

    py.xticks(np.arange(1, len(dataX)+1), ticksX, horizontalalignment='left', rotation=30)
    leg = py.legend()
    leg.draggable()
    py.xlabel(axisX)
    py.ylabel('time (s)')

    # Set X axis tick labels as rungs
    #print zip(dataX, dataY)
  
    py.draw()
    py.show()
    
    return

def interpret_msmc(top_param = program_parameters(), scaling_method = "2N0", year = 1, ylog10scale = False ):
    if scaling_method == "generation":
        scaling_method = "years"
        year = 1

    #msmc_para = program_parameters()
    #top_param  = msmc_para
    ####### This two lines should be outside, after called
    
    replicates = top_param.replicates
    nsample    = top_param.nsample
    case       = top_param.case
      
    dir_name = "msmc-experiment/" + case + "Samples" +`nsample`
    ms_param = param.ms_param_of_case(case)
    
    ms_param.plot(ylog10scale = ylog10scale, timescale = scaling_method, year = year)
    
    mu = ms_param.t / ms_param.seqlen / float(4) / ms_param.scaling_N0 

    for ith_run in range(replicates):
        ms_out_file_prefix = ms_param.case + \
                             "Samples" +`nsample` + \
                             "msdata"

        outputFile_name = dir_name + "/" + ms_out_file_prefix + ".final.txt"
        if os.path.isfile( outputFile_name ):
            #print outputFile_name
            outputFile = open ( outputFile_name, "r")
            tmp_time , tmp_lambda = [] , []        
            skipline = outputFile.readline()
           
            for line in outputFile:
                tmp_time.append( float(line.split()[1]) )
                tmp_lambda.append( float(line.split()[3]) )
            
            #print tmp_time, tmp_lambda
            if scaling_method == "years":
                time = [ ti / mu * year for ti in tmp_time]
            elif scaling_method == "2N0":
                time = [ ti / mu / ( 2* float(ms_param.scaling_N0) ) for ti in tmp_time]
            elif scaling_method == "4N0":
                time = [ ti / mu / ( 4* float(ms_param.scaling_N0)) for ti in tmp_time]
                
            time[0] = time[1]/float(100)
            time.append(time[-1]*float(100))
            yaxis_scaler = float(10000)
            pop = [ 1 / lambda_i / (2*mu) / yaxis_scaler for lambda_i in tmp_lambda]
            pop.insert(0, pop[0])              
            pylab.step(time, pop, color = "purple", linewidth=2.0)
            #print pop

            outputFile.close()
    pylab.savefig( dir_name + ".pdf" )            
    pylab.close()            

def graphical_output(account_balance_giro, account_balance_extra, depot_balance_extra, giro_data,
                     daily_account_balance):
    years = mdates.YearLocator()  # every year
    months = mdates.MonthLocator()
    weeks = mdates.WeekdayLocator()
    days = mdates.DayLocator()

    dateFormat = mdates.DateFormatter('%Y-%m-%d')
    yearFormat = mdates.DateFormatter('%Y-%m')

    fig = pl.figure()

    pl.xlabel("date")
    pl.ylabel("balance")
    ax1 = pl.subplot(121)
    pl.xlabel("date")
    pl.ylabel("account balance")
    pl.title("account balances")

    ax1.step(daily_account_balance[:, 0], daily_account_balance[:, 1])
    ax1.step(daily_account_balance[:, 0], daily_account_balance[:, 2])
    ax1.step(daily_account_balance[:, 0], daily_account_balance[:, 3])
    ax1.step(daily_account_balance[:, 0], daily_account_balance[:, 4])
    ax1.xaxis.set_major_locator(months)
    ax1.xaxis.set_major_formatter(yearFormat)
    # ax1.xaxis.set_minor_locator(weeks)
    ax1.format_xdata = mdates.DateFormatter('%Y-%m-%d')
    ax1.xaxis_date()
    pl.legend(["Giro Account", "Extra Account", "Depot", "Total"], loc="best")

    ax2 = pl.subplot(122)
    pl.xlabel("date")
    pl.ylabel("bookings")
    pl.title("giro account bookings")
    ax2.bar(giro_data[:, 0], giro_data[:, 4])
    # ax2.bar(extra_data[:,0], extra_data[:,4]) too much stuff going on
    ax2.xaxis.set_major_locator(months)
    ax2.xaxis.set_major_formatter(yearFormat)
    ax2.xaxis.set_minor_locator(weeks)
    ax2.format_xdata = mdates.DateFormatter('%Y-%m-%d')
    ax2.xaxis_date()
    fig.autofmt_xdate()

    fig = pl.figure()

    daily_account_balance_diff = daily_account_balance[30:, :] - daily_account_balance[:-30, :]
    daily_account_balance_diff[:, 0] = daily_account_balance[30:, 0]
    # print daily_account_balance_diff
    pl.xlabel("date")
    pl.ylabel("balance")
    pl.step(daily_account_balance_diff[:, 0], daily_account_balance_diff[:, -1])
    pl.axhline(y=0, xmin=0, xmax=1)
    pl.legend(["Total 30 day difference"], loc="best")

    pl.show()

def my_cdf(data):
    data.sort()
    n = len(data)
    y = []
    for i in range(n):
	y.append((float(i)+1)/n)

    p.step(data, y)
    p.title('CDF')
    p.show()        
    pass

def my_ccdf(data):
    # in-place sort
    data.sort()
    n = len(data)
    y = []
    for i in range(n):
	y.append(1-(float(i)+1)/n)
    p.step(data, y)
    p.ylabel('ccdf(x)')
    p.xlabel('x')
    p.title('CCDF')
    p.show()        
    pass

def plotFluxE(fluxVec, energyVec=None, fnameOut='fluxE', figNum=0, label='fluxE'):
    if not energyVec:
        # assume 10 energy grp structure by default
        energyVec = np.array([1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7])
    fluxVec = np.append(np.array([0]), fluxVec)
    plt.figure(figNum)
    plt.step(energyVec, fluxVec / np.sum(fluxVec), linewidth='4', label=label)
    plt.ylim([5e-3, 1])
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel("Energy (eV)")
    plt.ylabel("Log Flux (Arbitrary Scaling)")
    plt.legend()
    plt.savefig(fnameOut)

def measure_flexure(sky):
    ''' Measure expected (589.3 nm) - measured emission line in nm'''
    ll, ss = sky['nm'], sky['ph_10m_nm']

    pl.figure()
    pl.step(ll, ss)

    pix = SG.argrelmax(ss, order=20)[0]
    skynms = chebval(pix, sky['coefficients'])
    for s in skynms: pl.axvline(s)

    ixmin = np.argmin(np.abs(skynms - 589.3))
    dnm = 589.3 - skynms[ixmin]

    return dnm

def diCal_figure (cum_change, tmrca, position, prefix):
    tmrca = [x*2 for x in tmrca] # convert tmrca from unit of 4Ne to 2Ne
    site, absorptionTime = extract_diCal_time ( prefix )
    maxtime = max( max(absorptionTime), max(tmrca) )
    
    fig = pl.figure(figsize=(24,7)) 
    pl.plot(site, absorptionTime, color = "green", linewidth=3.0)
    pl.step(cum_change, tmrca, color = "red", linewidth=5.0)
    pl.plot(position, [maxtime*1.05]*len(position), "bo")    
    pl.axis([min(site), max(site) , 0, maxtime*1.1])
    pl.title("Dical Absorption time (Green)")
    pl.xlabel("Sequence base")
    pl.ylabel("TMRCA (2N0)")
    pl.savefig( prefix + "diCal_absorption_time" + ".png")
    pl.close()

def xsPlot(xs, energyVec=None, micro=True, label='xs1', style='-', fnameOut='xsPlt', figNum=4):
    if not energyVec:
        # assume 10 energy grp structure by default
        energyVec = np.array([1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7])
    xsVec = np.append(np.array([0]), xs)
    plt.figure(figNum)
    plt.step(energyVec, xsVec, linewidth='4', label=label, linestyle=style)
    plt.yscale('log')
    plt.xscale('log')
    plt.xlabel("Energy (eV)")
    if micro:
        plt.ylabel("Micro Xsection [barns]")
    else:
        plt.ylabel("Macro Xsection [1/cm]")
    plt.legend()
    plt.savefig(fnameOut)

def plot_attenuation(rmArray, f, label):
	N = rmArray.shape[0]
	l2 = (0.3 / f) ** 2
	factors = []
	for i, RM in enumerate(rmArray):
		atten = 0
		for j in range(i + 1, N, 1):
			atten += np.cos(2 * (rmArray[j] - RM) * l2)

		factor = (float(N) + (2. * atten)) / pow(float(N), 2.)
		factors.append(factor)

	Pv, v = np.histogram(factors, bins=10, density=True)
	print('Epsilon=', np.mean(factors), '+/-', np.std(factors))
	v = 0.5 * (v[1:] + v[:-1])
	pylab.step(v, Pv, label=label)

	return [np.mean(factors), np.std(factors)]