Python pylab.loglog函数代码示例

def logskew(d,floors = [],col='',divvypk=[]):
    maxd = max(d.flatten())
    if len(floors) == 0:
        floors = 1./maxd * 2.**N.arange(-5,5.01,1.)
    print floors

    for fl in floors:
        dcopy = 1.*d
        wltf = N.where(dcopy < fl)
        dcopy[wltf] = 0.*dcopy[wltf] + fl

        logrho = N.log(dcopy)
        var  = N.var(decreaseres(logrho,f=16).flatten())
        #skew = SS.skew(logrho.flatten())
        #print fl,'log:',var,'lin:',N.var(dcopy.flatten())
        print fl,'log:',var#,SS.skew(dcopy.flatten())

        k,pk = power.pk(logrho)

        if (len(divvypk) > 0):
            M.loglog(k,pk/divvypk/var,col)

    if (len(divvypk) == 0):
        return pk,pk/var
    else:
        return

def extrapolate(n, y, tolerance=1e-15, plot=False, call_show=True):
    "Extrapolate functional value Y from sequence of values (n, y)."

    # Make sure we have NumPy arrays
    n = array(n)
    y = array(y)

    # Create initial "bound"
    Y0 = 0.99*y[-1]
    Y1 = 1.01*y[-1]

    # Compute initial interior points
    phi = (sqrt(5.0) + 1.0) / 2.0
    Y2 = Y1 - (Y1 - Y0) / phi
    Y3 = Y0 + (Y1 - Y0) / phi

    # Compute initial values
    F0, e, nn, ee, yy = _eval(n, y, Y0)
    F1, e, nn, ee, yy = _eval(n, y, Y1)
    F2, e, nn, ee, yy = _eval(n, y, Y2)
    F3, e, nn, ee, yy = _eval(n, y, Y3)

    # Solve using direct search (golden ratio fraction)
    while Y1 - Y0 > tolerance:

        if F2 < F3:
            Y1, F1 = Y3, F3
            Y3, F3 = Y2, F2
            Y2 = Y1 - (Y1 - Y0) / phi
            F2, e, nn, ee, yy = _eval(n, y, Y2)
        else:
            Y0, F0 = Y2, F2
            Y2, F2 = Y3, F3
            Y3 = Y0 + (Y1 - Y0) / phi
            F3, e, nn, ee, yy = _eval(n, y, Y3)

        print Y0, Y1

    # Compute reference value
    Y = 0.5*(Y0 + Y1)

    # Print results
    print
    print "Reference value:", Y

    # Plot result
    if plot:
        pylab.figure()
        pylab.subplot(2, 1, 1)
        pylab.title("Reference value: %g" % Y)
        pylab.semilogx(n, y, 'b-o')
        pylab.semilogx(nn, yy, 'g--')
        pylab.subplot(2, 1, 2)
        pylab.loglog(n, e, 'b-o')
        pylab.loglog(nn, ee, 'g--')
        pylab.grid(True)
        if call_show:
            pylab.show()

    return Y

    def romanzuniga07(wavelength, AKs, makePlot=False):
        filters = ['J', 'H', 'Ks', '[3.6]', '[4.5]', '[5.8]', '[8.0]']
        wave =      np.array([1.240, 1.664, 2.164, 3.545, 4.442, 5.675, 7.760])
        A_AKs =     np.array([2.299, 1.550, 1.000, 0.618, 0.525, 0.462, 0.455])
        A_AKs_err = np.array([0.530, 0.080, 0.000, 0.077, 0.063, 0.055, 0.059])
        
        # Interpolate over the curve
        spline_interp = interpolate.splrep(wave, A_AKs, k=3, s=0)

        A_AKs_at_wave = interpolate.splev(wavelength, spline_interp)
        A_at_wave = AKs * A_AKs_at_wave

        if makePlot:
            py.clf()
            py.errorbar(wave, A_AKs, yerr=A_AKs_err, fmt='bo', 
                        markerfacecolor='none', markeredgecolor='blue',
                        markeredgewidth=2)

            # Make an interpolated curve.
            wavePlot = np.arange(wave.min(), wave.max(), 0.1)
            extPlot = interpolate.splev(wavePlot, spline_interp)
            py.loglog(wavePlot, extPlot, 'k-')

            # Plot a marker for the computed value.
            py.plot(wavelength, A_AKs_at_wave, 'rs',
                    markerfacecolor='none', markeredgecolor='red',
                    markeredgewidth=2)
            py.xlabel('Wavelength (microns)')
            py.ylabel('Extinction (magnitudes)')
            py.title('Roman Zuniga et al. 2007')


        return A_at_wave

    def plotppf(self,x=None,xmin=None,alpha=None,dolog=True,**kwargs):
        """
        Plots the power-law-predicted value on the Y-axis against the real
        values along the X-axis.  Can be used as a diagnostic of the fit 
        quality.
        """
        if not(xmin): xmin=self._xmin
        if not(alpha): alpha=self._alpha
        if not(x): x=numpy.sort(self.data[self.data>xmin])
        else: x=numpy.sort(x[x>xmin])

        # N = M^(-alpha+1)
        # M = N^(1/(-alpha+1))
        
        m0 = min(x)
        N = (1.0+numpy.arange(len(x)))[::-1]
        xmodel = m0 * N**(1/(1-alpha)) / max(N)**(1/(1-alpha))
        
        if dolog:
            pylab.loglog(x,xmodel,'.',**kwargs)
            pylab.gca().set_xlim(min(x),max(x))
            pylab.gca().set_ylim(min(x),max(x))
        else:
            pylab.plot(x,xmodel,'.',**kwargs)
        pylab.plot([min(x),max(x)],[min(x),max(x)],'k--')
        pylab.xlabel("Real Value")
        pylab.ylabel("Power-Law Model Value")

def plotData(filename):
    gatetime, allanvar = np.loadtxt(filename, comments="#", delimiter=",", unpack=True)
    pl.loglog(gatetime, allanvar, '.-')
    pl.xlabel("Gate time (s)")
    pl.ylabel("Allan deviation (Hz)")
    pl.xlim([min(gatetime), max(gatetime)])
    pl.show()

def plot_calibrated(time_len=128, chan=1):
    """ Plot calibrated as a function of time """
    nar.load_data() if plot_old else nar.take_data(time_len)
        
    cal   = (nar.ts_x_on + nar.ts_x_off) * Td / (nar.ts_x_on/nar.ts_x_off - 1)
    p_tot = (nar.ts_x_on + nar.ts_x_off)
    time_series  = cal[:, chan]
    p_avg        = np.average(p_tot[:, chan])
    time_series  = time_series / p_avg
    spec_series = np.abs(np.fft.rfft(time_series))
    
    #uncal = nar.ts_x_on
    #uncal = uncal[:, chan]
    #spec_uncal = np.abs(np.fft.rfft(uncal))
    
    t = np.cumsum(np.ones([time_len]))
    t = t /np.max(t) * total_time
    tu = np.cumsum(np.ones([time_len/2+1]))/total_time
    
    plt.subplot(211)
    #plt.plot(t, uncal, c=c[1])
    plt.plot(t, time_series, c=c[0])
    plt.xlabel("Time (s)")
    plt.ylabel("Calibrated signal (K)")
    plt.subplot(212)
    #plt.loglog(tu, spec_uncal, c=c[1])
    plt.loglog(tu, spec_series, c=c[0])

    plt.xlabel("Spectrum (Hz)")
    plt.ylabel("")
    plt.show()

def snr_mat_f(mchvec, reds, lum_dist, fmin, fmax, fvec, finteg, tobs, sn_f):
	''''''
	mch_fmat=np.transpose(np.tile(mchvec, (len(reds), len(fvec), finteg, 1) ), axes=(0,3,1,2))
	z_fmat=np.transpose(np.tile(reds, (len(mchvec), len(fvec), finteg, 1) ),axes=(3,0,1,2))
	f_fmat=np.transpose(np.tile(fvec, (len(reds), len(mchvec), finteg, 1) ), axes=(0,1,3,2))
	finteg_fmat=np.transpose(np.tile(np.arange(finteg), (len(reds), len(mchvec), len(fvec), 1) ), axes=(0,1,2,3))
	stshape=np.shape(z_fmat) #Standard shape of all matrices that I will use.
	DL_fmat=np.transpose(np.tile(lum_dist, (len(mchvec), len(fvec), finteg, 1) ),axes=(3,0,1,2)) #Luminosity distance in Mpc.
	flim_fmat=A8.f_cut(1./4., 2.*mch_fmat*2.**(1./5.))*1./(1.+z_fmat) #The symmetric mass ratio is 1/4, since I assume equal masses.
	flim_det=np.maximum(np.minimum(fmax, flim_fmat), fmin) #The isco frequency limited to the detector window.
	tlim_fmat=CM.tafter(mch_fmat, f_fmat, flim_fmat, z_fmat)
	#By construction, f_mat cannot be smaller than fmin or larger than fmax (which are the limits imposed by the detector).
	fmin_fmat=np.minimum(f_fmat, flim_det) #I impose that the minimum frequency cannot be larger than the fisco.
	fmaxobs_fmat=flim_det.copy()
	#fmaxobs_fmat=fmin_fmat.copy()
	fmaxobs_fmat[tobs<tlim_fmat]=CM.fafter(mch_fmat[tobs<tlim_fmat], z_fmat[tobs<tlim_fmat], f_fmat[tobs<tlim_fmat], tobs)
	fmax_fmat=np.minimum(fmaxobs_fmat, flim_det) #The maximum frequency (after an observation tobs) cannot exceed fisco or the maximum frequency of the detector.
	integconst=(np.log10(fmax_fmat)-np.log10(fmin_fmat))*1./(finteg-1)
	finteg_fmat=fmin_fmat*10**(integconst*finteg_fmat)
	sn_vec=sn_f(fvec)##########
	sn_fmat=sn_f(finteg_fmat) #Noise spectral density.
	#htilde_fmat=A8.htilde_f(1./4., 2.*mch_fmat*2**(1./5.), z_fmat, DL_fmat, f_fmat)
	htilde_fmat=A8.htilde_f(1./4., 2.*mch_fmat*2**(1./5.), z_fmat, DL_fmat, finteg_fmat)
	py.loglog(finteg_fmat[0,0,:,0],htilde_fmat[0,0,:,0]**2.)
	py.loglog(finteg_fmat[0,0,:,0],sn_fmat[0,0,:,0])
	snrsq_int_fmat=4.*htilde_fmat**2./sn_fmat #Integrand of the S/N square.
	snrsq_int_m_fmat=0.5*(snrsq_int_fmat[:,:,:,1:]+snrsq_int_fmat[:,:,:,:-1]) #Integrand at the arithmetic mean of the infinitesimal intervals.
	df_fmat=np.diff(finteg_fmat, axis=3) #Infinitesimal intervals.
	snr_full_fmat=np.sqrt(np.sum(snrsq_int_m_fmat*df_fmat,axis=3)) #S/N as a function of redshift, mass and frequency.
	fopt=fvec[np.argmax(snr_full_fmat, axis=2)] #Frequency at which the S/N is maximum, for each pixel of redshift and mass.
	snr_opt=np.amax(snr_full_fmat, axis=2) #Maximum S/N at each pixel of redshift and mass.
	snr_min=snr_full_fmat[:,:,0]
	return snr_opt

def nishiyama09(wavelength, AKs, makePlot=False):
    # Data pulled from Nishiyama et al. 2009, Table 1

    filters = ['V', 'J', 'H', 'Ks', '[3.6]', '[4.5]', '[5.8]', '[8.0]']
    wave =      np.array([0.551, 1.25, 1.63, 2.14, 3.545, 4.442, 5.675, 7.760])
    A_AKs =     np.array([16.13, 3.02, 1.73, 1.00, 0.500, 0.390, 0.360, 0.430])
    A_AKs_err = np.array([0.04,  0.04, 0.03, 0.00, 0.010, 0.010, 0.010, 0.010])

    # Interpolate over the curve
    spline_interp = interpolate.splrep(wave, A_AKs, k=3, s=0)

    A_AKs_at_wave = interpolate.splev(wavelength, spline_interp)
    A_at_wave = AKs * A_AKs_at_wave

    if makePlot:
        py.clf()
        py.errorbar(wave, A_AKs, yerr=A_AKs_err, fmt='bo', 
                    markerfacecolor='none', markeredgecolor='blue',
                    markeredgewidth=2)
        
        # Make an interpolated curve.
        wavePlot = np.arange(wave.min(), wave.max(), 0.1)
        extPlot = interpolate.splev(wavePlot, spline_interp)
        py.loglog(wavePlot, extPlot, 'k-')

        # Plot a marker for the computed value.
        py.plot(wavelength, A_AKs_at_wave, 'rs',
                markerfacecolor='none', markeredgecolor='red',
                markeredgewidth=2)
        py.xlabel('Wavelength (microns)')
        py.ylabel('Extinction (magnitudes)')
        py.title('Nishiyama et al. 2009')

    
    return A_at_wave

def ConvIndicator(X, Y, pct=0.1, fs=14, eqaxis=False):
    """
    Convergence indicator icon
    ==========================
    """
    if len(X)<2: raise Exception('at least 2 points are required')
    xx, yy   = log10(X), log10(Y)
    p        = polyfit(xx, yy, 1)
    m        = round(p[0])
    xx0, xx1 = min(xx), max(xx)
    yy0, yy1 = min(yy), max(yy)
    dxx, dyy = xx1-xx0, yy1-yy0
    xxm, yym = (xx0+xx1)/2.0, (yy0+yy1)/2.0
    xxl, xxr = xxm-pct*dxx, xxm+pct*dxx
    shift    = 0.5*pct*dxx*m
    xm,  ym  = 10.0**xxm, 10.0**(yym-shift)
    xl,  xr  = 10.0**xxl, 10.0**xxr
    yl,  yr  = 10.0**(yym+m*(xxl-xxm)-shift),10.0**(yym+m*(xxr-xxm)-shift)
    loglog(X, Y)
    #plot(xm, ym, 'ro')
    #plot(xl, yl, 'go')
    #plot(xr, yr, 'mo')
    points = array([[xl,yl],[xr,yl],[xr,yr]])
    gca().add_patch(Polygon(points, ec='k', fc='None'))
    xxR = xxm+1.2*pct*dxx
    xR  = 10.0**xxR
    text(xR, ym, '%g'%m, ha='left', va='center', fontsize=fs)
    if eqaxis: axis('equal')
    return m

def plot_point(rank1, style=None):
    if not style:
        style = "r--"
    # plot lines to point (rank1,pop1)
    pop1 = pop[rank1]
    P.loglog([rank1, rank1], [0.1, pop1], style)
    P.loglog([1, rank1], [pop1, pop1], style)

def main():
    plt.ion()

    fil = FletcherFilter()
    Niter = 12
    logp = plt.zeros((Niter,2))
    for k in range(Niter):
        while True:
            #print k
            p = plt.rand(2)
            if not fil.dominated(p):
                break
        logp[k] = p
        fil.add(p, 0.0, 0.0)
        ff = fil.values[fil.valid]
        ff = plt.r_[[[1e-6,1]], ff[plt.argsort(ff[:,0])], [[1,1e-6]]]
        ww = plt.zeros((ff.shape[0] * 2 - 1, 2))
        ww[::2] = ff
        ww[1::2,0] = ff[1:,0]
        ww[1::2,1] = ff[:-1,1]
        plt.loglog(ww[:,0], ww[:,1], '-')
    plt.loglog(logp[:,0], logp[:,1], 'ys-', lw=2)
    plt.axis([0,1,0,1])
    plt.axis('equal')
    plt.grid()
        
    code.interact()

def test_integrators(tmax=100,
                     Nsteps=1000,
                     x0=[0.0,1.0],
                     v0=[1.0,0.0]):
    mag = lambda x: numpy.sqrt(numpy.dot(x,x))
    fprime = lambda x: -x/mag(x)**3
    energy = lambda x,v: 0.5*(v**2).sum(1) - 1./numpy.sqrt((x**2).sum(1))

    pylab.figure(1)
    
    for (method,N) in ((Euler,Nsteps),
                       (DriftKick,Nsteps),
                       (LeapFrog,Nsteps),
                       (VelocityVerlet,Nsteps),
                       (RungeKutta,Nsteps/4)):
        t = numpy.linspace(0,tmax,N+1)
        x_t,v_t = method(fprime,x0,v0,tmax,N)
        E = energy(x_t,v_t)

        delta_E = abs((E-E[0])/E[0])
        
        pylab.loglog(t[1:],delta_E[1:],
                     label=method.__name__)
    
    pylab.ylim(1E-8,1E4)
    pylab.legend(loc=2)
    pylab.xlabel(r'$\mathdefault{t}$')
    pylab.ylabel(r'$\mathdefault{|\Delta E/E|}$')
    
    pylab.show()

def plot_convergence_data_from_file( labels ):
    import pickle
    l2_error = pickle.load( open( 'l2-error.p', 'rb' ) )
    inf_error = pickle.load( open( 'inf-error.p', 'rb' ) )
    num_pts = pickle.load( open( 'num-grid-points.p', 'rb' ) )


    max_error = 0.
    min_error = 1.
    for i in xrange( len( l2_error ) ):
        pylab.figure( 1 )
        pylab.loglog( num_pts[i], l2_error[i], '-o', label = labels[i] )


        pylab.figure( 2 )
        pylab.loglog( num_pts[i], inf_error[i], '-o', label = labels[i] )

        max_error = max( max_error, numpy.array( l2_error[i] ).max() )
        min_error = min( min_error, numpy.array( l2_error[i] ).min() )

    #matplotlib.rcParams.update({'font.size': 16})
    pylab.figure( 1 )
    pylab.legend()
    pylab.xlabel(r'Number of grid points')
    pylab.ylabel(r'$\lVert f - \hat{f}\rVert_{\ell_2}$')#,fontsize=16)
    pylab.ylim(min_error/5.,5.*max_error)
    pylab.savefig('genz-corner-peak-10d-5e-1c-quartic-decay-l2-convergence.eps',dpi=1200)

    pylab.figure( 2 )
    #pylab.title(r'$\int_{-1}^1 f(x)$', fontsize=10)
    pylab.legend()
    pylab.xlabel(r'Number of grid points')
    pylab.ylabel(r'$\lVert I[f] - Q[f]\rVert_{\ell_2}$')
    pylab.ylim(1e-16,5*max_error)

def expy(d):
    colors = ['b','g','r','c','m','y','k','b--','g--','r--','c--','m--','y--']
    multi = 2.**N.arange(2,3.1,1.)
    #multi = [77.]
    for i in range(len(multi)):
        dinflate = N.exp((d)*multi[i])-1.
        k,p = power.pk(dinflate)

        var = N.var((d)*multi[i])
        lognocor = ((N.exp(var)-1)*N.exp(var)/var)
        lognovar = (N.exp(var)-1)*N.exp(var)
        print 1+2*var,'lognocor:',lognocor
        #varlog = N.var(N.log(dinflate+1.).flatten())
        k,plog = power.pk(N.log(dinflate+1))
        print p[0]/plog[0], p[-1]/plog[-1]

        #M.subplot(121)
        M.loglog(k,p/plog/lognocor,colors[i])
        #M.subplot(122)

        #M.loglog(k,plog,colors[i])
        #M.loglog([k[0],k[-1]],[lognovar,lognovar],colors[i])
        #M.loglog(k,1./plog**mul,colors[i])
        

    
    #preal = M.load('mill/s63/pm.pnl.dat')
    #plogreal = M.load('mill/s63/plogm.pnl.dat')    
    #M.loglog(preal[:,0],preal[:,1]/plogreal[:,1],'b')

    #M.xlabel(r'$k\ [\rm{Mpc}/h]$',fontsize=20)
    #M.ylabel(r'$P_\delta(k)/P_{\log (1+\delta)}(k)$',fontsize=20)

    M.show()

def MisorientationScalingCollapseCompareInset(misses, bdlengths, labels, alpha=2.5):
    """ good values for alpha seem to be 4, but 2.5 for experiment """
    colors = ['b', 'r', 'g', 'y']

    pl.rcParams.update({'legend.fontsize': 14,
            'legend.columnspacing':1.2,
        })
    for i, (mis, label, bdlength) in enumerate(zip(misses, labels, bdlengths)):
        t = mis*bdlength/(mis*bdlength).mean()
        dx = 5./100.
        y,x = np.histogram(t, bins=np.linspace(0, 5, 100))
        x = (x[:-1]+x[1:])/2
        y = y.astype('float')/y.sum() / dx
        pl.plot(x, y, colors[i]+'o--', label=label)

        alpha = fit_alpha(x, y)
        xt = np.linspace(0,5, 1000)
        pl.plot(xt, scaling(alpha, xt), colors[i]+'-', label=r"Fit, $\alpha$ = %0.1f" % alpha)

    pl.xlabel(r"$\theta / \theta_{av}$")
    pl.ylabel(r"$\theta_{av}\,P(\theta, \theta_{av})$")
    pl.legend(loc='lower right')

    ax = pl.axes([0.52, 0.52, 0.35, 0.35])
    for i, (mis, label, bdlength) in enumerate(zip(misses, labels, bdlengths)):
        t = mis*bdlength/(mis*bdlength).mean()
        dx = 5./100.
        y,x = np.histogram(t, bins=np.linspace(0, 5, 100))
        x = (x[:-1]+x[1:])/2
        y = y.astype('float')/y.sum() / dx
        pl.plot(x, y, colors[i]+'o-', label=label)
    pl.loglog()

    return x, y