Python WLanalysis.smooth方法代码示例

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
def compute_GRF_PDF_ps_pk (r):
	'''for a convergence map with filename fn, compute the PDF and the power spectrum. sizedeg = 3.5**2, or 1.7**2'''
	print cosmo, r
	kmap = kmapGen(r)
	#kmap = load(CMBlensing_dir+'GRF_fidu/'+'GRF_fidu_%04dr.npy'%(r))
	
	i_arr = arange(len(sigmaP_arr))
	
	if not doGRF:
            kmap_smoothed = [WLanalysis.smooth(kmap, sigmaP) for sigmaP in sigmaP_arr]
            ps = WLanalysis.PowerSpectrum(kmap_smoothed[0])[1]

            PDF = [PDFGen(kmap_smoothed[i], PDFbin_arr[i]) for i in i_arr]
            peaks = [peaksGen(kmap_smoothed[i], peak_bins_arr[i]) for i in i_arr]

	###### generate GRF
	else:
            ps=0
            random.seed(r)
            GRF = (WLanalysis.GRF_Gen(kmap)).newGRF()
            #save(CMBlensing_dir+'GRF_fidu/'+'GRF_fidu_%04dr.npy'%(r), GRF)		
            #GRF = load(CMBlensing_dir+'GRF_fidu/'+'GRF_fidu_%04dr.npy'%(r))
            GRF_smoothed = [WLanalysis.smooth(GRF, sigmaP) for sigmaP in sigmaP_arr]
            PDF = [PDFGen(GRF_smoothed[i], PDFbin_arr[i]) for i in i_arr]
            peaks = [peaksGen(GRF_smoothed[i], peak_bins_arr[i]) for i in i_arr]
	#############

	return [ps,], PDF, peaks#, PDF_GRF, peaks_GRF

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
def iskew (i):
    print i
    ikmap_NL = kmapNL(i)
    ikmap_NOISY = kmapNOISY(i)
    skewness_NL = [skew(WLanalysis.smooth(ikmap_NL, ismooth).flatten() ) for ismooth in sigmaG_arr*PPA_NL] 
    skewness_NOISY = [skew(WLanalysis.smooth(ikmap_NOISY, ismooth).flatten() ) for ismooth in sigmaG_arr*PPA_NOISY]
    return [skewness_NL, skewness_NOISY]

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
def compute_GRF_PDF_ps_pk (cosmo, r, Gaus=0,sigmaG=8.0):
    kmap = FT2real(cosmo, r, Gaus=Gaus)
    ps = WLanalysis.PowerSpectrum(WLanalysis.smooth(kmap, 0.18), bins=bins)[1]#*2.0*pi/ell_arr**2
    if not filtered:
        kmap = WLanalysis.smooth(kmap, 2.93*sigmaG/8.0)
    PDF = PDFGen(kmap, PDFbins)
    peaks = peaksGen(kmap, peak_bins)
    return concatenate([ps, PDF, peaks])

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
	def plot_predict_maps_fcn(WxsigmaG):
		Wx, sigmaG = WxsigmaG
		mask0 = maskGen(Wx, 0.5, sigmaG)
		mask = WLanalysis.smooth(mask0, 5.0)
		
		kmap_lensing = kmap_lensing_Gen(Wx, sigmaG)
		kmap_predict = kmap_predict_Gen(Wx, sigmaG)
		bmode_lensing= bmode_lensing_Gen(Wx, sigmaG)
		
		mask_nan = mask0.copy()
		mask_nan[mask0==0]=nan
		
		#imshow(kmap_lensing*mask_nan, vmax=3*std(kmap_lensing), vmin=-2*std(kmap_lensing), origin = 'lower')
		#title('W%i kmap_lensing'%(Wx))
		#colorbar()
		#savefig(plot_dir+'kmap_W%i_sigmaG%s_lensing.jpg'%(Wx,sigmaG))
		#close()

		##imshow(kmap_predict, vmax=3*std(kmap_predict), vmin=-2*std(kmap_predict), origin = 'lower')
		imshow(kmap_predict*mask_nan, vmax=4*std(kmap_predict), vmin=0, origin = 'lower')
		#imshow(kmap_predict, origin = 'lower')
		title('W%i kmap_predict'%(Wx))
		colorbar()
		savefig(plot_dir+'kmap_L12_W%i_sigmaG%s_predict.jpg'%(Wx,sigmaG))
		close()

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
def KSmap(iinput):
	'''Input:
	i = ith zbin for zcut
	hl = 'hi' or 'lo' for higher/lower z of the zcut
	sigmaG: smoothing scale
	Wx = 1..4 of the field
	Output:
	smoothed KS map and galn map.
	'''
	Wx, sigmaG, i, hl = iinput
	print 'Wx, sigmaG, i, hl:', Wx, sigmaG, i, hl
	kmap_fn = cat_dir+'KS/W%i_KS_%s_%s_sigmaG%02d.fit'%(Wx, zbins[i],hl,sigmaG*10)
	galn_smooth_fn = cat_dir+'KS/W%i_galn_%s_%s_sigmaG%02d.fit'%(Wx, zbins[i],hl,sigmaG*10)
	
	isfile_kmap, kmap = WLanalysis.TestFitsComplete(kmap_fn, return_file = True)
	if isfile_kmap == False:
		Me1_fn = cat_dir+'Me_Mw_galn/W%i_Me1w_%s_%s.fit'%(Wx, zbins[i],hl)
		Me2_fn = cat_dir+'Me_Mw_galn/W%i_Me2w_%s_%s.fit'%(Wx, zbins[i],hl)
		Mw_fn = cat_dir+'Me_Mw_galn/W%i_Mwm_%s_%s.fit'%(Wx, zbins[i],hl)
		Me1 = WLanalysis.readFits(Me1_fn)
		Me2 = WLanalysis.readFits(Me2_fn)
		Mw = WLanalysis.readFits(Mw_fn)	
		Me1_smooth = WLanalysis.weighted_smooth(Me1, Mw, PPA=PPA512, sigmaG=sigmaG)
		Me2_smooth = WLanalysis.weighted_smooth(Me2, Mw, PPA=PPA512, sigmaG=sigmaG)
		kmap = WLanalysis.KSvw(Me1_smooth, Me2_smooth)
		WLanalysis.writeFits(kmap,kmap_fn)
	isfile_galn, galn_smooth = WLanalysis.TestFitsComplete(galn_smooth_fn, return_file = True)
	if isfile_galn == False:
		galn_fn = cat_dir+'Me_Mw_galn/W%i_galn_%s_%s.fit'%(Wx, zbins[i],hl)
		galn = WLanalysis.readFits(galn_fn)
		galn_smooth = WLanalysis.smooth(galn, sigma=sigmaG*PPA512)
		WLanalysis.writeFits(galn_smooth, galn_smooth_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
def iskew_GRF (i):    
    print i
    a=load('/Users/jia/weaklensing/CMBnonGaussian/colin_noisy/kappaMapTT_Gauss_10000sims/kappaMap%04dTT_3.pkl'%(i))
    areal = real(fftpack.ifft2(a))
    inorm = (2*pi*3.5/360.0)/(77.0**2)
    areal /= inorm 
    skewness_NOISY = [skew(WLanalysis.smooth(areal, ismooth).flatten() ) for ismooth in sigmaG_arr*PPA_NOISY]
    return skewness_NOISY

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
def maskGen_init (Wx, sigma_pix=10, JCH=0):
	galn = WLanalysis.smooth(load(cat_dir+'Me_Mw_galn/W%i_galn_zcut13.npy'%(Wx)),PPA512)
	if JCH:
		galn *= JCHPSmaskGen(Wx)
	else:
		galn *= PSmaskGen(Wx)## add point source mask for cmbl
	mask = zeros(shape=galn.shape)
	mask[10:-10,10:-10] = 1 ## remove edge 10 pixels
	idx = where(galn<0.5)
	mask[idx] = 0
	mask_smooth = WLanalysis.smooth(mask, sigma_pix)	
	######## print out fksy and fsky 2 ##########
	sizedeg = (sizes[Wx-1]/512.0)**2*12.0
	fsky = sum(mask_smooth)/sizes[Wx-1]**2*sizedeg/41253.0
	fsky2 = sum(mask_smooth**2)/sizes[Wx-1]**2*sizedeg/41253.0
	fmask = sum(mask_smooth)/sizes[Wx-1]**2
	fmask2 = sum(mask_smooth**2)/sizes[Wx-1]**2
	print 'W%i, fsky=%.8f, fsky2=%.8f, fmask=%.8f, fmask2=%.8f'%(Wx, fsky, fsky2, fmask,fmask2) 
	#############################################
	return mask_smooth#fsky, fsky2#

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
		def return_kappa_arr (Wx, sigmaG=sigmaG):
			mask = maskGen(Wx, 0.5, sigmaG)
			kmap_predict = kmap_predict_Gen(Wx, sigmaG)
			kmap_predict -= mean(kmap_predict)
			kmap_lensing = kmap_lensing_Gen(Wx, sigmaG)
			bmode = bmode_lensing_Gen(Wx, sigmaG)
			kproj_peak_mat = WLanalysis.peaks_mat(kmap_predict)
			#kproj_peak_mat = WLanalysis.peaks_mat(kmap_lensing)
			#kproj_peak_mat = WLanalysis.peaks_mat(bmode)
			idx_pos = (kproj_peak_mat!=0)&(~isnan(kproj_peak_mat))&(mask>0)
			kappa_proj = kmap_predict[idx_pos]
			kappa_lensing = kmap_lensing[idx_pos]
			kappa_bmode = bmode[idx_pos]
			
			######## do an overlay of peaks on top of convergence #######
			if sigmaG == 8.9:
				kmap_predict2 = kmap_predict_Gen(Wx, 5.3)
				mask2 = maskGen(Wx, 0.5, 5.3)
				kproj_peak_mat = WLanalysis.peaks_mat(kmap_predict2)
				kproj_peak_mat[mask2==0] = nan
				kproj_peak_mat[isnan(kproj_peak_mat)]=0
				peaksmooth = WLanalysis.smooth(kproj_peak_mat,10)
				kstd=std(kmap_lensing)
				#pstd=std(kmap_predict2)
				
				kmap_lensing[peaksmooth>2*std(peaksmooth)]=nan
				kmap_lensing[mask2==0]=-99
				kmap_predict2[peaksmooth>2*std(peaksmooth)]=nan
				f2=figure(figsize=(20,12))
				axx=f2.add_subplot(121)
				axy=f2.add_subplot(122)
				axx.imshow(kmap_lensing,origin='lower',vmin=-2*kstd,vmax=3*kstd,interpolation='nearest')
				#f2.colorbar()
				axx.set_title('lensing')
				axy.imshow(kmap_predict2,origin='lower',interpolation='nearest')
				#plt.colorbar(cax=axy)
				axy.set_title('predict')
				savefig(plot_dir+'peaks_location_W%s.jpg'%(Wx))
				close()
				
			return kappa_proj, kappa_lensing, kappa_bmode

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
		icat = icat[:,idx_k]
		
		#####################################
		f_Wx = WLanalysis.gnom_fun(center)	
		y, x = array(f_Wx(icat[:2]))
		weight = icat[3]
		#k = np.load(obsPK_dir+'kappa_predict_W%i.npy'%(Wx))#kappa_predict_Mmax2e15_W%i.npy
		A, galn = WLanalysis.coords2grid(x, y, array([k*weight, weight, k]), size=isize)
		Mkw, Mw, Mk = A
		###########################################
		
		for sigmaG in  (0.5, 1.0, 1.8, 3.5, 5.3, 8.9):
			print Wx, sigmaG
			
			mask0 = maskGen(Wx, 0.5, sigmaG)
			mask = WLanalysis.smooth(mask0, 5.0)
			################ make maps ######################
			kmap_predict = WLanalysis.weighted_smooth(Mkw, Mw, PPA=PPA512, sigmaG=sigmaG)
			kmap_predict*=mask
			np.save(obsPK_dir+'maps/kmap_W%i_predict_sigmaG%02d.npy'%(Wx, sigmaG*10), kmap_predict)
			###########################################
if plot_predict_maps:
	def plot_predict_maps_fcn(WxsigmaG):
		Wx, sigmaG = WxsigmaG
		mask0 = maskGen(Wx, 0.5, sigmaG)
		mask = WLanalysis.smooth(mask0, 5.0)
		
		kmap_lensing = kmap_lensing_Gen(Wx, sigmaG)
		kmap_predict = kmap_predict_Gen(Wx, sigmaG)
		bmode_lensing= bmode_lensing_Gen(Wx, sigmaG)
		

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
def k3Gen(r):
    print r
    kmap=kmapGen(r)
    kmap_smoothed=WLanalysis.smooth(kmap,sigmaG)
    k3=mean(kmap_smoothed**3)
    return k3

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
    cat_w0 = load('/Users/jia/weaklensing/CFHTLenS_downloads/CFHTLens_2015-08-18T04-37-45.npy').T
    RA, DEC, weight, Z_B,  MAG_i, MAG_y, MAG_u, MAG_g, MAG_r, MAG_z = cat_w0
    MAGI = amin(array([abs(MAG_y), abs(MAG_i)]),axis=0)    

    for Wx in range(1,5):
        print Wx
        center = centers[Wx-1]
        mask=xmaskGen(Wx)
        for cut in (22,23,24):
            #
            idx_Wx = where((RA<RAs[Wx-1][1])&(RA>RAs[Wx-1][0])&(MAGI>18)&(MAGI<cut)&(weight>0))[0]
            igaln = coords2grid_counts(array([RA,DEC]).T[idx_Wx], Wx)
            igaln=igaln/weightGen(Wx)
            igaln=igaln/mean(igaln[mask>0])-1
            igaln[mask<1]=0
            igaln_smooth=WLanalysis.smooth(igaln,1.0)
            save('/Users/jia/weaklensing/multiplicative/referee/galn_W%i_cut%i_LensfitWNonzero.npy'%(Wx,cut),igaln_smooth)

if create_dndz_LensfitWNonzero:
    z_arr = arange(.025,3.5,0.05)
    Pz = load('/Users/jia/weaklensing/CFHTLenS_downloads/2015-09-290split/Pz_xac.npy')
    #Pz = (np.core.defchararray.replace(Pz,',','')).astype(float)
    RA, DEC, weight, ZB = load('/Users/jia/weaklensing/CFHTLenS_downloads/2015-09-290split/ra_dec_weightz_xac.npy').T
    cat_w0 = load('/Users/jia/weaklensing/CFHTLenS_downloads/CFHTLens_2015-08-18T04-37-45.npy').T
    master_ra, master_dec, master_weight, Z_B,  MAG_i, MAG_y, MAG_u, MAG_g, MAG_r, MAG_z = cat_w0
    MAGI = amin(array([abs(MAG_y), abs(MAG_i)]),axis=0)
    
    from pylab import *
    figure()

    for cut in (22,23,24):

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
        #savefig(CMBNG_dir+'plot_official/plot_noiseless_%s_diff.pdf'%(['PDF','peaks'][j-1]))
        #close()

if plot_sample_noiseless_noisy_map:
    import pickle
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    #kmap_noiseless=WLanalysis.readFits(CMBNG_dir+'test_maps/WLconv_z1100.00_0001r.fits')
    #kmap_noiseless_8arcmin=WLanalysis.smooth(kmap_noiseless,78.01904762)
    kmap_noiseless_8arcmin=load(CMBNG_dir+'test_maps/kmap_noiseless_8arcmin_1r.npy')
    #### get noisy map ########
    FTmap = pickle.load(open(CMBNG_dir+'colin_noisy/kappaMapTT_10000sims/kappaMap0000TT_3.pkl'))
    areal = real(fftpack.ifft2(FTmap))
    inorm = (2*pi*3.5/360.0)/(77.0**2)
    areal /= inorm
    kmap_noisy=areal
    kmap_noisy_8arcmin = WLanalysis.smooth(kmap_noisy, 2.93)
    
    f=figure(figsize=(12,5))
    
    i=1
    for kmap in [kmap_noiseless_8arcmin, kmap_noisy_8arcmin]:
        ax=subplot(1,2,i)
        #imshow(kmap_noiseless_8arcmin,vmin=-0.1,vmax=0.1,extent=(0,3.5,0,3.5))
        im=ax.imshow(kmap,vmin=-3*std(kmap_noiseless_8arcmin),vmax=3*std(kmap_noiseless_8arcmin),extent=(0,3.5,0,3.5),cmap='PuOr')
        ax.annotate(r"$\rm{%s}$"%(['noiseless','noisy'][i-1]), 
                    xy=(0.05, 0.9), xycoords='axes fraction',fontsize=24,
                    bbox={'facecolor':'thistle', 'alpha':0.5})
        ax.set_xlabel(r"$\rm{deg}$",fontsize=22)
        ax.set_ylabel(r"$\rm{deg}$",fontsize=22)
        ax.tick_params(labelsize=16)
        divider=make_axes_locatable(ax)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import smooth [as 别名]
	for fn in os.listdir(tSZ_dir+'planck/'):
		if fn[-3:]=='txt':
			print fn
			print fn[:-3]+'npy'
			full_fn = tSZ_dir+'planck/'+fn
			imap=WLanalysis.txt2map_fcn(full_fn, offset = False)
		if fn[-3:]=='npy':
			imap = load(tSZ_dir+'planck/'+fn)
		if 'mask' in fn:
			imshow(imap, origin='lower', vmin=0, vmax=1)
		elif '857' in fn:
			imshow(imap, origin='lower')
		elif 'JCH_ymap' in fn:
			imshow(imap, origin='lower', vmin=-2e-5, vmax=2e-5)
		elif 'GARY' in fn:
			imap = WLanalysis.smooth(imap, 2.5*4)
			imshow(imap, origin='lower', vmin=-2e-5, vmax=2e-5)
		else:
			#imshow(imap, origin='lower', vmin=-3*std(imap), vmax=3*std(imap))
			imshow(imap, origin='lower', vmin=-3e-6, vmax=3e-6)
		title(fn[:-4])
		colorbar()
		savefig(plot_dir+fn[:-3]+'jpg')
		close()

########### cross correlation	
cat_dir = '/Users/jia/weaklensing/CFHTLenS/catalogue/'

kmapGen = lambda i: load(cat_dir+'kmap_W%i_sigma10_zcut13.npy'%(i))

PSmaskGen = lambda i: np.load(tSZ_dir+'planck/PSmaskHFIall_flipper2048_CFHTLS_W%i.npy'%(i))