Python FramedPlot.xrange方法代码示例

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def plot_boss_geometry(color=None, colorwheel=None, plt=None, width=1, show=True,
                       region=None):
    """
    Plot the boundaries in the boss_survey.par file
    """
    import esutil as eu
    import biggles
    from biggles import FramedPlot, Curve

    bg = read_boss_geometry()

    if plt is None:
        plt = FramedPlot()
        plt.xlabel=r'$\lambda$'
        plt.ylabel=r'$\eta$'

    if color is not None:
        colors = [color]*len(bg)
    elif colorwheel is not None:
        colors = colorwheel
    else:
        colors = ['red','blue','green','magenta','navyblue','seagreen',
                  'firebrick','cadetblue','green4']
        

    for i in xrange(len(bg)):
        b = bg[i]
        color = colors[i % len(colors)]
        c = eu.plotting.bbox( b['clambdaMin'], b['clambdaMax'], b['cetaMin'], b['cetaMax'],
                             color=color, width=width)
        plt.add(c)

    if region == 'ngc':
        plt.yrange = [-40.,50.]
        plt.xrange = [-80.,80.]
    elif region == 'sgc':
        plt.yrange = [105.,165.]
        plt.xrange = [-60.,60.]
    else:
        plt.yrange = [-40.,165.]
        plt.xrange = [-80.,80.]


    plt.aspect_ratio = (plt.yrange[1]-plt.yrange[0])/(plt.xrange[1]-plt.xrange[0])

    if show:
        plt.show()

    return plt

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def test_fit_nfw_dsig(rmin=0.01):

    from biggles import FramedPlot,Points,SymmetricErrorBarsY,Curve
    omega_m=0.25
    z=0.25
    n = lensing.nfw.NFW(omega_m, z)

    r200 = 1.0
    c = 5.0

    rmax = 5.0
    log_rmin = log10(rmin)
    log_rmax = log10(rmax)
    npts = 25
    logr = numpy.linspace(log_rmin,log_rmax,npts)
    r = 10.0**logr

    ds = n.dsig(r, r200, c)
    # 10% errors
    dserr = 0.1*ds
    ds += dserr*numpy.random.standard_normal(ds.size)

    guess = numpy.array([r200,c],dtype='f8')
    # add 10% error to the guess
    guess += 0.1*guess*numpy.random.standard_normal(guess.size)

    res = fit_nfw_dsig(omega_m, z, r, ds, dserr, guess)

    r200_fit = res['r200']
    r200_err = res['r200_err']

    c_fit = res['c']
    c_err = res['c_err']

    print 'Truth:'
    print '    r200: %f' % r200
    print '       c: %f' % c
    print 'r200_fit: %f +/- %f' % (r200_fit,r200_err)
    print '   c_fit: %f +/- %f' % (c_fit,c_err)
    print 'Cov:'
    print res['cov']

    logr = numpy.linspace(log_rmin,log_rmax,1000)
    rlots = 10.0**logr
    yfit = n.dsig(rlots,r200_fit,c_fit)

    plt=FramedPlot()
    plt.add(Points(r,ds,type='filled circle'))
    plt.add(SymmetricErrorBarsY(r,ds,dserr))
    plt.add(Curve(rlots,yfit,color='blue'))

    plt.xlabel = r'$r$ [$h^{-1}$ Mpc]'
    plt.ylabel = r'$\Delta\Sigma ~ [M_{sun} pc^{-2}]$'

    plt.xrange = [0.5*rmin, 1.5*rmax]
    plt.yrange = [0.5*(ds-dserr).min(), 1.5*(ds+dserr).max()]

    plt.xlog=True
    plt.ylog=True
    plt.show()

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def test_interp_hybrid():
    """
    Send y0,y1 with one or both less than zero to test the
    hybrid offset scheme
    """
    slope = -2.0
    #xvals = 0.1+linspace(0.0,8.0,9)
    xvals = 10.0**linspace(0.0,1.0,10)
    yvals = xvals**slope
    #yerr = 0.5*yvals
    #yerr = sqrt(yvals)
    yerr = yvals.copy()
    yerr[:] = 0.05

    #xfine = 0.1+linspace(0.0,8.0,1000)
    xfine = 10.0**linspace(0.0,1.0,1000)
    yfine = xfine**slope 

    #yerr = yvals.copy()
    #yerr[:] = 2

    plt=FramedPlot()
    plt.xrange = [0.5*xvals.min(), 1.5*xvals.max()]
    plt.xlog=True
    #plt.ylog=True
    plt.add(Points(xvals,yvals,type='filled circle',size=1))
    plt.add(Curve(xfine,yfine,color='blue'))
    #w=where1( (yvals-yerr) > 1.e-5 )
    #plt.add(SymmetricErrorBarsY(xvals[w],yvals[w],yerr[w]))
    plt.add(SymmetricErrorBarsY(xvals,yvals,yerr))

    # make points in between consecutive xvals,yvals so we
    # can use hybrid 2-point function
    xi = numpy.zeros(xvals.size-1,dtype='f8')
    yi = numpy.zeros(xi.size,dtype='f8')
    for i in xrange(xi.size):

        logx = (log10(xvals[i+1])+log10(xvals[i]))/2.0
        xi[i] = 10.0**logx
        yi[i],amp,slope,off = interp_hybrid(xvals[i], xvals[i+1], 
                                            yvals[i], yvals[i+1],
                                            yerr[i], yerr[i+1], 
                                            xi[i],more=True)
        
        print 'amp:',amp
        print 'slope:',slope
        print 'off:',off

    print xvals
    print xi
    print yi
    plt.add( Points(xi, yi, type='filled circle', size=1, color='red'))

    plt.show()

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def plot_sub_pixel(ellip,theta, show=False):
    import biggles
    from biggles import PlotLabel,FramedPlot,Table,Curve,PlotKey,Points
    from pcolors import rainbow

    f=subpixel_file(ellip,theta,'fits')
    data = eu.io.read(f)
    colors = rainbow(data.size,'hex')

    pltSigma = FramedPlot()
    pltSigma.ylog=1
    pltSigma.xlog=1

    curves=[]
    for j in xrange(data.size):
        sigest2 = (data['Irr'][j,:] + data['Icc'][j,:])/2

        pdiff = sigest2/data['sigma'][j]**2 -1
        nsub=numpy.array(data['nsub'][j,:])

        #pc = biggles.Curve(nsub, pdiff, color=colors[j])
        pp = Points(data['nsub'][j,:], pdiff, type='filled circle',color=colors[j])

        pp.label = r'$\sigma: %0.2f$' % data['sigma'][j]
        curves.append(pp)
        pltSigma.add(pp)
        #pltSigma.add(pc)
        #pltSigma.yrange=[0.8,1.8]
        #pltSigma.add(pp)


    c5 = Curve(linspace(1,8, 20), .005+zeros(20))
    pltSigma.add(c5)

    key=PlotKey(0.95,0.95,curves,halign='right',fontsize=1.7)
    key.key_vsep=1

    pltSigma.add(key)
    pltSigma.xlabel='N_{sub}'
    pltSigma.ylabel=r'$\sigma_{est}^2  /\sigma_{True}^2 - 1$'

    lab=PlotLabel(0.05,0.07,r'$\epsilon: %0.2f \theta: %0.2f$' % (ellip,theta),halign='left')
    pltSigma.add(lab)

    pltSigma.yrange = [1.e-5,0.1]
    pltSigma.xrange = [0.8,20]
    if show:
        pltSigma.show()

    epsfile=subpixel_file(ellip,theta,'eps')

    print("Writing eps file:",epsfile)
    pltSigma.write_eps(epsfile)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
    def make_plot(self):

        bindata=self.bindata
        bin_field=self.bin_field
        plt=FramedPlot()
        
        plt.uniform_limits=1
        plt.xlog=True
        #plt.xrange=[0.5*bindata[bin_field].min(), 1.5*bindata[bin_field].max()]
        plt.xrange=self.s2n_range
        plt.xlabel=bin_field
        plt.ylabel = r'$\gamma$'
        if self.yrange is not None:
            plt.yrange=self.yrange

        xdata=bindata[bin_field]
        xerr=bindata[bin_field+'_err']

        if self.shnum in sh1exp:
            g1exp=zeros(xdata.size)+sh1exp[self.shnum]
            g2exp=zeros(xdata.size)+sh2exp[self.shnum]
            g1exp_plt=Curve(xdata, g1exp)
            g2exp_plt=Curve(xdata, g2exp)
            plt.add(g1exp_plt)
            plt.add(g2exp_plt)


        xerrpts1 = SymmetricErrorBarsX(xdata, bindata['g1'], xerr)
        xerrpts2 = SymmetricErrorBarsX(xdata, bindata['g2'], xerr)

        type='filled circle'
        g1color='blue'
        g2color='red'
        g1pts = Points(xdata, bindata['g1'], type=type, color=g1color)
        g1errpts = SymmetricErrorBarsY(xdata, bindata['g1'], bindata['g1_err'], color=g1color)
        g2pts = Points(xdata, bindata['g2'], type=type, color=g2color)
        g2errpts = SymmetricErrorBarsY(xdata, bindata['g2'], bindata['g2_err'], color=g2color)

        g1pts.label=r'$\gamma_1$'
        g2pts.label=r'$\gamma_2$'

        key=biggles.PlotKey(0.9,0.5,[g1pts,g2pts],halign='right')

        plt.add( xerrpts1, g1pts, g1errpts )
        plt.add( xerrpts2, g2pts, g2errpts )
        plt.add(key)

        labels=self.get_labels()

        plt.add(*labels)
        plt.aspect_ratio=1

        self.plt=plt

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def plot_dsig_one(r, dsig, dsigerr, **kw):
    """
    plot delta sigma

    useful if adding to an existing plot

    parameters
    ----------
    r: array
        radius
    dsig: array
        delta sigma
    dsigerr: array
        error on delta sigma
    """

    from biggles import FramedPlot

    nbin=1
    _set_biggles_defs(nbin)

    visible=kw.get('visible',True)
    xlog=kw.get('xlog',True)
    ylog=kw.get('ylog',True)
    aspect_ratio=kw.get('aspect_ratio',1)

    is_ortho=kw.get('is_ortho',False)

    plt=kw.get('plt',None)

    if plt is None:
        plt=FramedPlot()
        plt.aspect_ratio=aspect_ratio
        plt.xlog=xlog
        plt.ylog=ylog

        plt.xlabel = LABELS['rproj']
        if is_ortho:
            plt.ylabel = LABELS['osig']
        else:
            plt.ylabel = LABELS['dsig']

    xrng, yrng = _add_dsig_to_plot(plt, r, dsig, dsigerr, **kw)
    plt.xrange=xrng
    plt.yrange=yrng

    if visible:
        plt.show()

    return plt

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
    def make_frac_plot(self):
        if self.shnum not in sh1exp:
            raise ValueError("you must know the expected value")

        if sh1exp[self.shnum] != 0:
            gfield='g1'
            gtrue=sh1exp[self.shnum]
        elif sh2exp[self.shnum] != 0:
            gfield='g2'
            gtrue=sh2exp[self.shnum]
        else:
            raise ValueError("all expected are listed as zero")

        bindata=self.bindata
        bin_field=self.bin_field
        plt=FramedPlot()
        
        plt.title=self.get_title()

        plt.xlog=True
        plt.xrange=[0.5*bindata[bin_field].min(), 1.5*bindata[bin_field].max()]
        plt.xlabel=bin_field
        ylabel=r'$\Delta \gamma/\gamma$'
        plt.ylabel = ylabel

        xdata=bindata[bin_field]

        zero=zeros(xdata.size)
        zero_plt=Curve(xdata, zero)
        plt.add(zero_plt)

        xfill=[xdata.min(), xdata.max()]

        plt.add( biggles.FillBetween(xfill, [0.004,0.004], 
                                     xfill, [-0.004,-0.004],
                                     color='grey80'))


        gfrac = bindata[gfield]/gtrue-1
        gfrac_err = bindata[gfield+'_err']/gtrue
        type='filled circle'
        color='blue'
        gpts = Points(xdata, gfrac, type=type, color=color)
        gerrpts = SymmetricErrorBarsY(xdata, gfrac, gfrac_err,color=color)

        plt.add( gpts, gerrpts )

        self.plt=plt

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def plot_drho(comb=None, r=None, drho=None, drhoerr=None, 
              color='black',type='filled circle',
              nolabel=False, noshow=False, minval=1.e-5,
              aspect_ratio=1):
    """
    This one stands alone. 
    """

    if comb is not None:
        r=comb['rdrho']
        drho=comb['drho']
        drhoerr=comb['drhoerr']
    else:
        if r is None or drho is None or drhoerr is None:
            raise ValueError("Send a combined struct or r,drho,drhoerr")


    plt=FramedPlot()
    plt.aspect_ratio=aspect_ratio
    plt.xlog=True
    plt.ylog=True

    if not nolabel:
        plt.xlabel = r'$r$ [$h^{-1}$ Mpc]'
        plt.ylabel = r'$\delta\rho ~ [M_{\odot} pc^{-3}]$'


    od=add_to_log_plot(plt, r, drho, drhoerr, 
                       color=color, 
                       type=type,
                       minval=minval)

    # for drho we need even broader yrange
    plt.xrange = od['xrange']

    yr=od['yrange']
    plt.yrange = [0.5*yr[0], 3*yr[1]]

    if not noshow:
        plt.show()
    od['plt'] = plt
    return od

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def plot_mass(comb=None, r=None, mass=None, masserr=None, 
              color='black',type='filled circle',
              nolabel=False, noshow=False, minval=1.e11,
              aspect_ratio=1):

    if comb is not None:
        r=comb['rmass']
        mass=comb['mass']
        masserr=comb['masserr']
    else:
        if r is None or mass is None or masserr is None:
            raise ValueError("Send a combined struct or r,mass,masserr")


    plt=FramedPlot()
    plt.aspect_ratio=aspect_ratio
    plt.xlog=True
    plt.ylog=True

    if not nolabel:
        plt.xlabel = r'$r$ [$h^{-1}$ Mpc]'
        plt.ylabel = r'$M(<r) ~ [h^{-1} M_{\odot}]$'


    od=add_to_log_plot(plt, r, mass, masserr, 
                       color=color, 
                       type=type,
                       minval=minval)

    plt.xrange = od['xrange']
    plt.yrange = od['yrange']

    if not noshow:
        plt.show()
    od['plt'] = plt
    return od

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
    def plot_coverage(self, region='both', show=True, dops=True):
        import biggles
        from biggles import FramedPlot, Points, PlotKey


        l = self.read()
        w,=numpy.where( (l['ra'] >= 0.0) & (l['ra'] <= 360.0) )
        if w.size != l.size:
            print("threw out",l.size-w.size,"with bad ra")
            l=l[w]

        llam,leta = eu.coords.eq2sdss(l['ra'],l['dec'])
        maskflags = l['maskflags']

        lammin,lammax = (-70.,70.)

        if region=='ngc':
            # a bit high to make room for the legend
            emin,emax=(-40,60)

            biggles.configure('screen','width', 1800)
            biggles.configure('screen','height', 1140)

        elif region=='sgc':
            emin,emax=(100,165)
            biggles.configure('screen','width', 1800)
            biggles.configure('screen','height', 1140)
        else:
            emin,emax=(-40,165)
            biggles.configure('screen','width', 1140)
            biggles.configure('screen','height', 1140)

        wl=where1((leta > emin) & (leta < emax))
        llam=llam[wl]
        leta=leta[wl]
        maskflags=maskflags[wl]


        plt=FramedPlot()
        plt.xlabel=r'$\lambda$'
        plt.ylabel=r'$\eta$'

        print("adding all lenses")

        type = 'filled circle'
        symsize=0.2

        allp = Points(llam,leta,type=type,size=symsize)
        allp.label='all'
        plt.add(allp)

        wquad = es_sdsspy.stomp_maps.quad_check(maskflags)
        print("adding quad pass")
        quadp = Points(llam[wquad],leta[wquad],type=type,color='red',size=symsize)
        quadp.label = 'quad good'
        plt.add(quadp)

        fakepoints = eu.plotting.fake_filled_circles(['all','quad good'],['black','red'])
        key=PlotKey(0.95,0.95,fakepoints,halign='right')
        plt.add(key)


        es_sdsspy.stomp_maps.plot_boss_geometry(color='blue',plt=plt,show=False)

        xrng = (lammin, lammax)
        yrng = (emin, emax)
        plt.xrange = xrng
        plt.yrange = yrng
        plt.aspect_ratio = (yrng[1]-yrng[0])/float(xrng[1]-xrng[0])


        if show:
            plt.show()

        if dops:
            d = lensing.files.sample_dir(type='lcat',sample=self['sample'])
            d = os.path.join(d,'plots')
            if not os.path.exists(d):
                os.makedirs(d)
            epsfile = os.path.join(d, 'lcat-%s-%s-coverage.eps' % (self['sample'],region) )
            print("Writing to eps file:",epsfile)
            plt.write_eps(epsfile)
        return plt

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def doplot(run, st, s2n_field, setname, s2n_range, sratio_range, psfnums, nobj):
    tab=Table(2,1)

    color1='blue'
    color2='red'

    m1pts=Points(st['s2n'],st['m1'],type='filled circle',color=color1)
    m1errpts=SymmetricErrorBarsY(st['s2n'],st['m1'],st['m1_err'],color=color1)
    m2pts=Points(st['s2n'],st['m2'],type='filled circle',color=color2)
    m2errpts=SymmetricErrorBarsY(st['s2n'],st['m2'],st['m2_err'],color=color2)

    c1pts=Points(st['s2n'],st['c1'],type='filled circle',color=color1)
    c1errpts=SymmetricErrorBarsY(st['s2n'],st['c1'],st['c1_err'],color=color1)
    c2pts=Points(st['s2n'],st['c2'],type='filled circle',color=color2)
    c2errpts=SymmetricErrorBarsY(st['s2n'],st['c2'],st['c2_err'],color=color2)

    m1pts.label=r'$\gamma_1$'
    m2pts.label=r'$\gamma_2$'

    key=PlotKey(0.9,0.2,[m1pts,m2pts], halign='right')

    xrng=array( [0.5*s2n_range[0], 1.5*s2n_range[1]])
    cyrng=get_symmetric_range(st['c1'],st['c1_err'],st['c2'],st['c2_err'])
    myrng=get_symmetric_range(st['m1'],st['m1_err'],st['m2'],st['m2_err'])

    mplt=FramedPlot()
    cplt=FramedPlot()
    mplt.xlog=True
    cplt.xlog=True

    mplt.xrange=xrng
    cplt.xrange=xrng
    #mplt.yrange=myrng
    #cplt.yrange=cyrng
    mplt.yrange=[-0.15,0.15]
    cplt.yrange=[-0.01,0.01]

    mplt.xlabel=s2n_field
    mplt.ylabel='m'
    cplt.xlabel=s2n_field
    cplt.ylabel='c'

    zplt=Curve(xrng, [0,0])
    mallow=Curve(xrng, [-0.004, 0.004])
    callow=Curve(xrng, [-0.0004, 0.0004])

    mallow=FillBetween(xrng, [0.004,0.004], 
                       xrng, [-0.004,-0.004],
                       color='grey80')
    callow=FillBetween(xrng, [0.0004,0.0004], 
                       xrng, [-0.0004,-0.0004],
                       color='grey80')



    labels=get_labels(run,
                      psfnums,
                      setname,
                      nobj,
                      sratio_range)

    mplt.add( mallow, zplt, m1pts, m1errpts, m2pts, m2errpts, key )
    mplt.add(*labels)

    cplt.add( callow, zplt, c1pts, c1errpts, c2pts, c2errpts )

    tab[0,0] = mplt
    tab[1,0] = cplt

    tab.show()

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
    def plot_coverage_bybin(self, binner, region='both', 
                            show=True, dops=True, rand=None):
        import pcolors
        import biggles
        import converter
        from biggles import FramedPlot, Points, PlotKey


        orig = self.read_original()
        lcat = self.read()

        all_clam,all_ceta = eu.coords.eq2sdss(orig['ra'],orig['dec'])

        l = orig[lcat['zindex']]
        clam,ceta = eu.coords.eq2sdss(lcat['ra'],lcat['dec'])

        clammin,clammax = (-70.,120.)

        if region=='ngc':
            # a bit high to make room for the legend
            emin,emax=(-40,60)

            biggles.configure('screen','width', 1800)
            biggles.configure('screen','height', 1140)

            clammin,clammax = (-70.,120.)

        elif region=='sgc':
            emin,emax=(105,165)
            biggles.configure('screen','width', 1800)
            biggles.configure('screen','height', 1140)
            clammin,clammax = (-50.,90.)
        else:
            emin,emax=(-40,165)
            biggles.configure('screen','width', 1140)
            biggles.configure('screen','height', 1140)
            clammin,clammax = (-70.,120.)

        wl=where1((all_ceta > emin) & (all_ceta < emax))
        all_clam=all_clam[wl]
        all_ceta=all_ceta[wl]

        wl=where1((ceta > emin) & (ceta < emax))
        clam=clam[wl]
        ceta=ceta[wl]
        l=l[wl]


        plt=FramedPlot()
        plt.xlabel=r'$\lambda$'
        plt.ylabel=r'$\eta$'
        xrng = (clammin, clammax)
        yrng = (emin, emax)
        plt.xrange = xrng
        plt.yrange = yrng


        print("adding all lenses")

        type = 'filled circle'
        symsize=0.2
        colors = pcolors.rainbow(binner['nbin'],'hex')


        if rand is not None:
            clam_r,ceta_r = eu.coords.eq2sdss(rand['ra'],rand['dec'])
            wl=where1((ceta_r > emin) & (ceta_r < emax))
            clam_r=clam_r[wl]
            ceta_r=ceta_r[wl]
            rp = Points(clam_r, ceta_r, type='dot', size=0.2)
            plt.add(rp)
        
        size_min=0.2
        size_max=4

        sizes=[]
        minlambda = l['lambda_zred'].min()
        maxlambda = l['lambda_zred'].max()
        for i in xrange(binner['nbin']):
            w=binner.select_bin(l, i)
            mlam=l['lambda_zred'][w].mean()
            # scale 0 to 1
            sz=(mlam-minlambda)/maxlambda
            # now scale size
            sz = size_min + sz*(size_max-size_min)
            sizes.append(sz)

        all_plots=[]
        labels=[]
        #for i in xrange(binner['nbin']):
        for i in reversed(xrange(binner['nbin'])):
            w=binner.select_bin(l, i)

            #points = Points(clam[w], ceta[w],type=type,size=symsize, color=colors[i])
            points = Points(clam[w], ceta[w],type=type,size=sizes[i], color=colors[i])
            labels.append(binner.bin_label(i))

            plt.add(points)


#.........这里部分代码省略.........

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def plot_nfw_lin_fits_byrun(run, name, npts=100, prompt=False, 
                            withlin=True,
                            ymin=0.01, ymax=2000.0):
    """

    This should be made not specific for the m-z splits we
    used on the sims

    """
    conf = lensing.files.cascade_config(run)
    if withlin:
        ex='lin'
        nex='lin'
    else:
        nex=''
        ex=''
    d = lensing.sample_read(type='fit',sample=run, name=name, extra=ex)

    omega_m = conf['cosmo_config']['omega_m']

    rravel = d['r'].ravel()
    xrange = [0.5*rravel.min(), 1.5*rravel.max()]

    #for i in xrange(d.size):
    i=0
    for dd in d:

        zrange = dd['z_range']
        mrange = dd['m200_range']

        if dd['rrange'][0] > 0:
            log_rmin = log10(dd['rrange'][0])
            log_rmax = log10(dd['rrange'][1])
        else:
            log_rmin = log10(dd['r'][0])
            log_rmax = log10(dd['r'][-1])
        rvals = 10.0**linspace(log_rmin,log_rmax,npts)

        plt = FramedPlot()  
        lensing.plotting.add_to_log_plot(plt, dd['r'],dd['dsig'],dd['dsigerr'])

        z = dd['z_mean']
        fitter = lensing.fit.NFWBiasFitter(omega_m,z,rvals,withlin=withlin)

        if withlin:
            yfit = fitter.nfw_lin_dsig(rvals, dd['r200_fit'],dd['c_fit'],dd['B_fit'])
            yfit_nfw = fitter.nfw.dsig(rvals,dd['r200_fit'],dd['c_fit'])
            yfit_lin = fitter.lin_dsig(rvals,dd['B_fit'])

            yfit = where(yfit < 1.e-5, 1.e-5, yfit)
            yfit_lin = where(yfit_lin < 1.e-5, 1.e-5, yfit_lin)

            cyfit = Curve(rvals,yfit,color='blue')
            cyfit_nfw = Curve(rvals,yfit_nfw,color='red')
            cyfit_lin = Curve(rvals,yfit_lin,color='orange')

            cyfit.label = 'Best Fit'
            cyfit_nfw.label = 'NFW'
            cyfit_lin.label = 'linear'

            key=PlotKey(0.1,0.3,[cyfit,cyfit_nfw,cyfit_lin])
            plt.add(cyfit,cyfit_nfw,cyfit_lin,key)
        else:
            yfit_nfw = fitter.nfw.dsig(rvals,dd['r200_fit'],dd['c_fit'])
            cyfit_nfw = Curve(rvals,yfit_nfw,color='blue')
            plt.add(cyfit_nfw)

        zlab='%0.2f < z < %0.2f' % (zrange[0],zrange[1])
        plt.add(PlotLabel(0.7,0.8,zlab))
        ll = (log10(mrange[0]),log10(mrange[1]))
        mlab = r'$%0.2f < logM_{200} < %0.2f$' % ll
        plt.add(PlotLabel(0.7,0.9,mlab))

        #yrange = [ymin,(dd['dsig']+dd['dsigerr']).max()*1.5]
        yrange = [ymin,ymax]
        plt.xrange = xrange
        plt.yrange = yrange
        plt.xlog=True
        plt.ylog=True
        plt.xlabel = r'$r$ [$h^{-1}$ Mpc]'
        plt.ylabel = r'$\Delta\Sigma ~ [M_{sun} pc^{-2}]$'
        plt.aspect_ratio=1
        if prompt:
            plt.show()
            rinput = raw_input('hit a key: ')
            if rinput == 'q':
                return
        else:
            d = lensing.files.lensbin_plot_dir(run,name)
            if not os.path.exists(d):
                os.makedirs(d)
            epsfile=path_join(d,'desmocks-nfw%s-fit-%02i.eps' % (nex,i))
            print 'Writing epsfile:',epsfile
            plt.write_eps(epsfile)
        i += 1

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]

#.........这里部分代码省略.........
        else:
            print('  doing meane')
            be1 = eu.stat.Binner(field_data, self['e1'][w], weights=weights)
            be2 = eu.stat.Binner(field_data, self['e2'][w], weights=weights)
            ylabel = r'$<e>$'


        # regular hist for display
        print("  regular fixed binsize hist")
        xm,xe,xstd=eu.stat.wmom(field_data, weights, sdev=True)
        #hxmin = xm-4.0*xstd
        #hxmax = xm+4.0*xstd
        bsize = xstd/5.
        hist = eu.stat.histogram(field_data, binsize=bsize, 
                                 weights=weights, more=True)


        print("  hist  e1, nperbin: ",nperbin)
        be1.dohist(nperbin=nperbin, min=fmin, max=fmax)
        #be1.dohist(nbin=nbin, min=fmin, max=fmax)
        print("  stats e1")
        be1.calc_stats()
        print("  hist  e2, nperbin: ",nperbin)
        be2.dohist(nperbin=nperbin, min=fmin, max=fmax)
        #be2.dohist(nbin=nbin, min=fmin, max=fmax)
        print("  stats e2")
        be2.calc_stats()



        plt = FramedPlot()

        if xrng is not None:
            plt.xrange=xrng
        else:
            if field == 'R':
                plt.xrange=[0.29,1.01]

        if yrng is not None:
            plt.yrange=yrng
            ymin = yrng[0]
            ymax = 0.8*yrng[1]
        else:
            ymin = min( be1['wymean'].min(),be2['wymean'].min() )
            ymax = 0.8*max( be1['wymean'].max(),be2['wymean'].max() )

        # this is a histogram-like object
        ph = eu.plotting.make_hist_curve(hist['low'], hist['high'], hist['whist'], 
                                         ymin=ymin, ymax=ymax, color='grey50')
        plt.add(ph)

        p1 = Points( be1['wxmean'], be1['wymean'], 
                    type='filled circle', color='blue')
        p1err = SymErrY( be1['wxmean'], be1['wymean'], be1['wyerr2'], color='blue')
        p1.label = r'$e_1$'

        p2 = Points( be2['wxmean'], be2['wymean'], 
                    type='filled circle', color='red')
        p2.label = r'$e_2$'
        p2err = SymErrY( be2['wxmean'], be2['wymean'], be2['wyerr2'], color='red')

        key = PlotKey(0.1,0.9, [p1,p2])
        plt.add(p1, p1err, p2, p2err, key)

        if self.camcol != 'any' and field == 'R' and plot_type=='meane':
            order=3

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import xrange [as 别名]
def test_fit_nfw_lin_dsig(rmin=0.01):

    from biggles import FramedPlot,Points,SymmetricErrorBarsY,Curve,PlotKey
    omega_m=0.25
    z=0.25


    r200 = 1.0
    c = 5.0
    B=10.0

    rmax = 50.0
    log_rmin = log10(rmin)
    log_rmax = log10(rmax)
    npts = 30
    r = 10.0**linspace(log_rmin,log_rmax,npts)

    fitter = NFWBiasFitter(omega_m,z,r)
    ds = fitter.dsig(r, r200, c, B)
    # 10% errors
    dserr = 0.1*ds
    ds += dserr*numpy.random.standard_normal(ds.size)

    guess = numpy.array([r200,c,B],dtype='f8')
    # add 10% error to the guess
    guess += 0.1*guess*numpy.random.standard_normal(guess.size)

    res = fitter.fit(ds,dserr,guess, more=True)

    r200_fit=res['r200']
    r200_err = res['r200_err']
    c_fit=res['c']
    c_err = res['c_err']
    B_fit=res['B']
    B_err = res['B_err']


    print 'Truth:'
    print '    r200: %f' % r200
    print '       c: %f' % c
    print '       B: %f' % B
    print 'r200_fit: %f +/- %f' % (r200_fit,r200_err)
    print '   c_fit: %f +/- %f' % (c_fit,c_err)
    print '   B_fit: %f +/- %f' % (B_fit,B_err)
    print 'Cov:'
    print res['cov']

 
    rfine = 10.0**linspace(log_rmin,log_rmax,100)
    fitter2 = NFWBiasFitter(omega_m,z,rfine)

    yfit = fitter2.dsig(rfine, r200_fit, c_fit, B_fit)
    yfit_nfw = fitter2.nfw.dsig(rfine, r200_fit, c_fit)
    yfit_lin = fitter2.lin_dsig(rfine,B_fit)

    plt=FramedPlot()
    plt.add(Points(r,ds,type='filled circle'))
    plt.add(SymmetricErrorBarsY(r,ds,dserr))

    cyfit = Curve(rfine,yfit,color='blue')
    cyfit_nfw = Curve(rfine,yfit_nfw,color='red')
    cyfit_lin = Curve(rfine,yfit_lin,color='orange')

    cyfit.label = 'Best Fit'
    cyfit_nfw.label = 'NFW'
    cyfit_lin.label = 'linear'

    key=PlotKey(0.1,0.3,[cyfit,cyfit_nfw,cyfit_lin])
    plt.add(cyfit,cyfit_nfw,cyfit_lin,key)

    plt.xlabel = r'$r$ [$h^{-1}$ Mpc]'
    plt.ylabel = r'$\Delta\Sigma ~ [M_{sun} pc^{-2}]$'

    plt.xrange = [0.5*rmin, 1.5*rmax]
    plt.yrange = [0.5*(ds-dserr).min(), 1.5*(ds+dserr).max()]

    plt.xlog=True
    plt.ylog=True
    plt.show()