Python FramedPlot.write_eps方法代码示例

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
    def plot_masserr(self, filename):
        import pcolors
        data=eu.io.read(filename)
        nr200 = data.size
        nb = data['B'][0].size

        plt = FramedPlot()
        colors = pcolors.rainbow(nr200, 'hex')
        clist = []
        for ri in xrange(nr200):
            r200 = data['r200'][ri]
            m200 = data['m200'][ri]

            p = Points(data['B'][ri], (m200-data['m200meas'][ri])/m200, 
                       type='filled circle', color=colors[ri])
            crv = Curve(data['B'][ri], (m200-data['m200meas'][ri])/m200, 
                        color=colors[ri])
            crv.label = 'r200: %0.2f' % r200
            clist.append(crv)
            plt.add(p,crv)

        key = PlotKey(0.1,0.4,clist)
        plt.add(key)
        plt.xlabel=r'$\Omega_m \sigma_8^2 D(z)^2 b(M,z)$'
        plt.ylabel = r'$(M_{200}^{true}-M)/M_{200}^{true}$'

        epsfile = filename.replace('.rec','.eps')
        print 'writing eps:',epsfile
        plt.write_eps(eu.ostools.expand_path(epsfile))

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
    def plot_radec(self, type):
        """
        ra/dec plot of all points and the matched points
        """
        import biggles
        import converter
        from biggles import FramedPlot,Points

        print()
        dir=self.plotdir()
        if not os.path.exists(dir):
            os.makedirs(dir)

        psfile = self.plotfile(type)

        orig = self.read_original(type)
        mat = self.read_matched(type)

        plt=FramedPlot() 

        symsize = 2
        if type == 'sdss' or type == 'other':
            symsize = 0.25

        plt.add( Points(orig['ra'],orig['dec'],type='dot', size=symsize) )

        plt.add( Points(mat['ra'],mat['dec'],type='dot',
                        color='red', size=symsize) )
        plt.xlabel = 'RA'
        plt.ylabel = 'DEC'
        print("Writing eps file:", psfile)
        plt.write_eps(psfile)
        converter.convert(psfile, dpi=120, verbose=True)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
def plot_nfwfits_byrun(run, name, prompt=False):
    conf = lensing.files.read_config(run)
    d = lensing.sample_read(type='fit', sample=run, name=name)
    omega_m = conf['omega_m']


    rvals = numpy.linspace(d['r'].min(), d['r'].max(),1000)
    for i in xrange(d.size):
        plt = FramedPlot()  
        lensing.plotting.add_to_log_plot(plt, 
                                          d['r'][i],
                                          d['dsig'][i],
                                          d['dsigerr'][i])

        z = d['z_mean'][i]
        n = lensing.nfw.NFW(omega_m, z)
        yfit = n.dsig(rvals, d['r200_fit'][i],d['c_fit'][i])
        plt.add(Curve(rvals,yfit,color='blue'))
        plt.xlog=True
        plt.ylog=True
        plt.xlabel = r'$r$ [$h^{-1}$ Mpc]'
        plt.ylabel = r'$\Delta\Sigma ~ [M_{sun} pc^{-2}]$'
        if prompt:
            plt.show()
            raw_input('hit a key: ')
        else:
            epsfile='/home/esheldon/tmp/plots/desmocks-nfwfit-%02i.eps' % i
            print 'Writing epsfile:',epsfile
            plt.write_eps(epsfile)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [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 write_eps [as 别名]
    def plot_ellip_vs_input(self, show=False):
        '''
        Plot the measured ellip as a function of the input for sigma_index=0
        which is a reasonably large object
        '''
        import biggles
        from biggles import PlotLabel,FramedPlot,Table,Curve,PlotKey,Points
        from pcolors import rainbow
        import pprint

        data = self.read()
        w=where1(data['sigma_index'] == 10)
        data = data[w]
        
        e1_input, e2_input, Tinput = util.mom2ellip(data['Irr_input'],
                                                    data['Irc_input'],
                                                    data['Icc_input'])
        e1_meas, e2_meas, Tinput = util.mom2ellip(data['Irr_meas'],
                                                  data['Irc_meas'],
                                                  data['Icc_meas'])

        einput = sqrt(e1_input**2 + e2_input**2)
        emeas = sqrt(e1_meas**2 + e2_meas**2)

        plt=FramedPlot()

        p = Points(einput,emeas, type='filled circle')
        plt.add(p)
        plt.xlabel=r'$\epsilon_{in}$'
        plt.ylabel=r'$\epsilon_{meas}$'

        sig=sqrt((data['Irr_meas'][0]+data['Icc_meas'][0])/2)
        lab1=PlotLabel(0.1,0.9,self.model, halign='left')
        lab2=PlotLabel(0.1,0.8,r'$\sigma: %0.2f$' % sig, halign='left')
        plt.add(lab1,lab2)


        einput.sort()
        c = Curve(einput, einput, color='red')
        c.label = r'$\epsilon_{input} = \epsilon_{meas}$'
        key=PlotKey(0.95,0.07,[c], halign='right')
        plt.add(c)
        plt.add(key)
        if show:
            plt.show()

        epsfile=self.epsfile('ellip-vs-input')
        print("Writing eps file:",epsfile)
        plt.write_eps(epsfile)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
    def plot_size_vs_input(self, show=False):
        '''
        Plot recovered size vs input for ellip=0, which is ellip_index=0
        '''

        import biggles
        from biggles import PlotLabel,FramedPlot,Table,Curve,PlotKey,Points
        from pcolors import rainbow
        import pprint

        data = self.read()
        w=where1(data['ellip_index'] == 0)
        data = data[w]

        siginput = sqrt(data['Irr_input'])
        sigmeas  = sqrt(data['Irr_meas'])


        pars=numpy.polyfit(siginput, sigmeas, 1)
        print("offset:",pars[1])
        print("slope: ",pars[0])
        print("IGNORING OFFSET")

        plt=FramedPlot()

        p = Points(siginput,sigmeas, type='filled circle')
        plt.add(p)
        plt.xlabel=r'$\sigma_{in}$'
        plt.ylabel=r'$\sigma_{meas}$'

        lab=PlotLabel(0.1,0.9,self.model)
        plt.add(lab)

        yfit2=pars[0]*siginput
        cfit2=Curve(siginput, yfit2, color='steel blue')
        cfit2.label = r'$%0.2f \sigma_{in}$' % pars[0]

        plt.add( cfit2 )

        key=PlotKey(0.95,0.07,[cfit2], halign='right')
        plt.add(key)
        if show:
            plt.show()

        epsfile=self.epsfile('size-vs-input')
        print("Writing eps file:",epsfile)
        plt.write_eps(epsfile)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
    def compare_all_other(self, type, show=True):
        
        fdict=self.all_other_fdict(type)

        # this is the original file.  It has the redshifts
        orig = zphot.weighting.read_training(fdict['origfile'])

        # this is the outputs
        num = zphot.weighting.read_num(fdict['numfile1'])

        # this is the weights file
        weights = zphot.weighting.read_training(fdict['wfile2'])

        # recoverable set
        w_recoverable = where1(num['num'] > 0)
        # this is actually the indexes back into the "orig" file
        w_keep = num['photoid'][w_recoverable]

        # get the z values for these validation objects
        zrec = orig['z'][w_keep]

        binsize=0.0314
        valid_dict = histogram(zrec, min=0, max=1.1, binsize=binsize, more=True)
        plt=FramedPlot()

        vhist = valid_dict['hist']/(float(valid_dict['hist'].sum()))
        pvhist=biggles.Histogram(vhist, x0=valid_dict['low'][0], binsize=binsize)
        pvhist.label = 'truth'

        weights_dict = histogram(weights['z'], min=0, max=1.1, binsize=binsize,
                                 weights=weights['weight'], more=True)
        whist = weights_dict['whist']/weights_dict['whist'].sum()
        pwhist=biggles.Histogram(whist, x0=weights_dict['low'][0], 
                                 binsize=binsize, color='red')
        pwhist.label = 'weighted train'

        key = PlotKey(0.6,0.6,[pvhist,pwhist])
        plt.add(pvhist,pwhist,key)

        plt.add( biggles.PlotLabel(.8, .9, type) )

        plt.write_eps(fdict['zhistfile'])
        converter.convert(fdict['zhistfile'],dpi=90,verbose=True)
        if show:
            plt.show()

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
    def plot_ellip_vs_field(self, field, rmag_max=21.8, fmin=None, fmax=None, nbin=20, nperbin=50000,
                            yrange=None, show=True):
        self.load_data()

        w=where1((self['cmodelmag_dered_r'] > 18.0) & (self['cmodelmag_dered_r'] < rmag_max) )

        if w.size == 0:
            print("no good objects")
            return

        weights = 1.0/(0.32**2 + self['uncer_rg'][w]**2)

        if field == 'psf_fwhm':
            field_data = self['psf_fwhm'][w]
            fstr = 'PSF FWHM (arcsec)'
        elif field == 'psf_sigma':
            field_data = self['psf_sigma'][w]
            fstr = r'$\sigma_{PSF}$'
        elif field == 'R_rg':
            field_data = self['r_rg'][w]
            fstr = 'R_rg'
        else:
            field_data = self[field][w]
            fstr=field

        print("Plotting mean e for field:",field)

        fstr = fstr.replace('_','\_')

        be1 = eu.stat.Binner(field_data, self['e1_rg'][w], weights=weights)
        be2 = eu.stat.Binner(field_data, self['e2_rg'][w], weights=weights)

        print("  hist  e1")
        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")
        be2.dohist(nperbin=nperbin, min=fmin, max=fmax)
        #be2.dohist(nbin=nbin, min=fmin, max=fmax)
        print("  stats e2")
        be2.calc_stats()

        plt = FramedPlot()
        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.8, 0.9, [p1,p2])
        plt.add(p1, p1err, p2, p2err, key)

        if field != 'cmodelmag_dered_r':
            rmag_lab = PlotLabel(0.1,0.05,'rmag < %0.2f' % rmag_max, halign='left')
            plt.add(rmag_lab)

        plab = PlotLabel(0.1,0.1, 'CH+RM', halign='left')
        plt.add(plab)

        plt.xlabel = r'$'+fstr+'$'
        plt.ylabel = r'$<e>$'

        if yrange is not None:
            plt.yrange=yrange

        if show:
            plt.show()
        epsfile = self.plotfile(field, rmag_max)
        print("  Writing eps file:",epsfile)
        plt.write_eps(epsfile)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [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 write_eps [as 别名]

#.........这里部分代码省略.........
            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)


        labels.reverse()
        fakepoints = eu.plotting.fake_filled_circles(labels, colors)
        key=PlotKey(0.95,0.95,fakepoints,halign='right',size=1.5)
        plt.add(key)

        plt.aspect_ratio = (yrng[1]-yrng[0])/float(xrng[1]-xrng[0])

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

        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-coverage-bybin.eps' % self['sample'])
            if rand is not None:
                epsfile=epsfile.replace('.eps','-withrand.eps')
            if region in ['sgc','ngc']:
                epsfile=epsfile.replace('.eps','-%s.eps' % region)

            print("Writing to eps file:",epsfile)
            plt.write_eps(epsfile)
            print("converting to png")
            converter.convert(epsfile, dpi=300)
        return plt

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

#.........这里部分代码省略.........
                                 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
            print("  getting poly order",order)
            coeff1 = numpy.polyfit(be1['wxmean'], be1['wymean'], order)
            poly1=numpy.poly1d(coeff1)
            coeff2 = numpy.polyfit(be2['wxmean'], be2['wymean'], order)
            poly2=numpy.poly1d(coeff2)

            ps1 = Curve( be1['wxmean'], poly1(be1['wxmean']), color='blue')
            ps2 = Curve( be2['wxmean'], poly2(be2['wxmean']), color='red')
            plt.add(ps1,ps2)

            polyf = self.R_polyfile(self.camcol, rmag_max)
            out={'coeff_e1':list([float(c) for c in coeff1]),
                 'coeff_e2':list([float(c) for c in coeff2])}
            print("    -- Writing poly coeffs to:",polyf)
            eu.io.write(polyf,out)
        if field != 'rmag':
            rmag_lab = \
                PlotLabel(0.1,0.05,'%0.2f < rmag < %0.2f' % (rmag_min,rmag_max), 
                          halign='left')
            plt.add(rmag_lab)

        procrun_lab = PlotLabel(0.1,0.1,
                            'procrun: %s filter: %s' % (self.procrun, self.band), 
                            halign='left')
        plt.add(procrun_lab)
        cy=0.9
        if self.run != 'any':
            run_lab = PlotLabel(0.9,0.9, 'run: %06i' % self.run, halign='right')
            plt.add(run_lab)
            cy=0.8
        if self.camcol != 'any':
            run_lab = PlotLabel(0.9,cy, 'camcol: %i' % self.camcol, halign='right')
            plt.add(run_lab)



        plt.xlabel = r'$'+fstr+'$'
        plt.ylabel = ylabel


        if show:
            plt.show()
        epsfile = self.plotfile(field, rmag_max, plot_type=plot_type)
        eu.ostools.makedirs_fromfile(epsfile, verbose=True)
        print("  Writing eps file:",epsfile)
        plt.write_eps(epsfile)

        converter.convert(epsfile, verbose=True)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
def test_dsig(run, ntest, catdict=None, overwrite=False):

    outdir='/home/users/esheldon/www/tmp/test-shear'
    if catdict is None:
        catdict = load_test_data(run)

    conf = lensing.files.read_config(run)
    c = eu.cosmology.Cosmo(omega_m=conf['omega_m'])

    lcat = catdict['l']
    lout = catdict['lout']
    scat = catdict['s']

    xl=catdict['xl']
    yl=catdict['yl']
    zl=catdict['zl']
    xs=catdict['xs']
    ys=catdict['ys']
    zs=catdict['zs']

    m200_min = 3.e14
    wl = where1(lcat['m200'] > m200_min)

    rmin = 0.1
    #rmax = 1.0
    rmax = 0.5
    for ii in xrange(ntest):
        i = wl[ii]

        f = path_join(outdir,'lens-%04i.eps' % i)
        if not os.path.exists(f) or overwrite:
            print 'lens:',i
            print 'getting lens angular diameter distance'
            zlens   = lcat['z'][i]
            ralens  = lcat['ra'][i]
            declens = lcat['dec'][i]
            DL    = c.Da(0.0, zlens)[0]
            print '    zl: ',zlens
            print '    DL: ',DL

            maxangle = rmax/DL
            minangle = rmin/DL
            cos_maxangle = cos(maxangle)
            cos_minangle = cos(minangle)

            cosdis = xl[i]*xs + yl[i]*ys + zl[i]*zs

            w=where1( cosdis > cos_maxangle )

            print 'Number within %0.2f arcmin: %s/%s' \
                    % (maxangle*180/numpy.pi*60,w.size,scat.size)

            if w.size > 0:
                w2=where1( cosdis[w] < cos_minangle )
                if w2.size > 0:
                    w=w[w2]

                    # radians
                    r = arccos(cosdis[w])

                    # get total gamma instead bothering to project to tangent.
                    # Let's just see how it looks
                    #gamma = sqrt( s['g1'][w]**2 + s['g2'][w]**2 )
                    
                    # both radians
                    rr,theta = gcirc(ralens,
                                     declens,
                                     scat['ra'][w],
                                     scat['dec'][w],
                                     getangle=True)
                    rdiff = r-rr
                    print 'maxx diff rad:',numpy.abs(rdiff).max()
                    print theta[0:10]

                    r*=DL
                    rr*=DL
                    cos2theta = cos(2*theta)
                    sin2theta = sin(2*theta)
                    gammaT = -(scat['g1'][w]*cos2theta + scat['g2'][w]*sin2theta);

                    zsrc = scat['z'][w]

                    plt = FramedPlot()
                    p = Points(zsrc,gammaT*r,type='filled circle')
                    plt.add(p)
                    plt.xlabel = 'source z'
                    plt.ylabel = 'gamma*r (Mpc)'
                    plt.add(PlotLabel(0.2,0.9,'zlens: %0.2f' % zlens))
                    plt.add(PlotLabel(0.2,0.8,'m200: %0.2e' % lcat['m200'][i]))

                    print 'Writing to:',f
                    plt.write_eps(f)

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [as 别名]
    def compare_same_same(self, type, show=True):
        """
        Use the id from the validation set to go back and get the
        z for those objects.  Then plot histograms for comparision.

        read in all file
        read in validation set
            take recoverable subset based on num file
        Get z info for these points from the all file

        plot the histgram of actual validation set redshifts
        overplot the histgram of weighted redshifts

        Then bin by true validation set redshift and plot the
            ztrue - <z>
        Where <z> is the expectation value of z based on the p(z)
            <z> = integral( z*p(z) )/integral( p(z) )
        That will be noisy
        """
        
        fdict=self.same_same_fdict(type)

        # this is the original file
        all = zphot.weighting.read_training(fdict['origfile'])

        # this is the validation set, for which the "photoid" field
        # is actually an id pointing back into "all"
        # we take version 1 and will demand num > 0
        valid = zphot.weighting.read_photo(fdict['photofile'])
        num = zphot.weighting.read_num(fdict['numfile1'])


        # this is the weights file
        weights = zphot.weighting.read_training(fdict['wfile2'])

        # recoverable set
        w_recoverable = where1(num['num'] > 0)
        # this is actually the indexes back into the "all" file
        w_keep = num['photoid'][w_recoverable]

        # get the z values for these validation objects
        zvalid = all['z'][w_keep]

        binsize=0.0314
        valid_dict = histogram(zvalid, min=0, max=1.1, binsize=binsize, more=True)
        plt=FramedPlot()

        vhist = valid_dict['hist']/(float(valid_dict['hist'].sum()))
        pvhist=biggles.Histogram(vhist, x0=valid_dict['low'][0], binsize=binsize)
        pvhist.label = 'validation'

        weights_dict = histogram(weights['z'], min=0, max=1.1, binsize=binsize,
                                 weights=weights['weight'], more=True)
        whist = weights_dict['whist']/weights_dict['whist'].sum()
        pwhist=biggles.Histogram(whist, x0=weights_dict['low'][0], 
                                 binsize=binsize, color='red')
        pwhist.label = 'weighted train'

        key = PlotKey(0.6,0.6,[pvhist,pwhist])
        plt.add(pvhist,pwhist,key)

        plt.add( biggles.PlotLabel(.8, .9, type) )

        plt.write_eps(fdict['zhistfile'])
        converter.convert(fdict['zhistfile'],dpi=90,verbose=True)
        if show:
            plt.show()

# 需要导入模块: from biggles import FramedPlot [as 别名]
# 或者: from biggles.FramedPlot import write_eps [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