Python WLanalysis.readFits方法代码示例

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def Bmode(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 'Bmode - Wx, sigmaG, i, hl:', Wx, sigmaG, i, hl
	bmap_fn = cat_dir+'KS/W%i_Bmode_%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, bmap = WLanalysis.TestFitsComplete(bmap_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)
		### Bmode conversion is equivalent to
		### gamma1 -> gamma1' = -gamma2
		### gamma2 -> gamma2' = gamma1
		bmap = WLanalysis.KSvw(-Me2_smooth, Me1_smooth)
		WLanalysis.writeFits(bmap,bmap_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def TestCrossCorrelate (Wx, zcut, sigmaG):
	'''Input: 
	Wx - one of the W1..W4 field (= 1..4) 
	zcut - redshift cut between KS background galaxies and forground cluster probe
	sigmaG - smoothing
	Output:
	ell_arr, CCK, CCB
	'''
	galn_hi = WLanalysis.readFits(test_dir+'W%i_galn_%s_hi_sigmaG%02d.fit'%(Wx,zcut,sigmaG*10))
	galn_lo = WLanalysis.readFits(test_dir+'W%i_galn_%s_lo_sigmaG%02d.fit'%(Wx,zcut,sigmaG*10))
	galn_cut = 0.5*0.164794921875 #5gal/arcmin^2*arcmin^2/pix, arcmin/pix = 12.0*60**2/512.0**2 = 
	bmap = WLanalysis.readFits(test_dir+'W%i_Bmode_%s_hi_sigmaG%02d.fit'%(Wx,zcut,sigmaG*10))
	kmap = WLanalysis.readFits(test_dir+'W%i_KS_%s_hi_sigmaG%02d.fit'%(Wx,zcut,sigmaG*10))
	mask = where(galn_hi<galn_cut)
	bmap[mask]=0
	kmap[mask]=0
	edges=linspace(5,100,11)
	ell_arr, CCB = WLanalysis.CrossCorrelate (bmap,galn_lo,edges=edges)
	ell_arr, CCK = WLanalysis.CrossCorrelate (kmap,galn_lo,edges=edges)
	f=figure(figsize=(8,6))
	ax=f.add_subplot(111)
	ax.plot(ell_arr, CCB, 'ro',label='B-mode')
	ax.plot(ell_arr, CCK, 'bo', label='KS')
	legend()
	#ax.set_xscale('log')
	ax.set_xlabel('ell')
	ax.set_ylabel(r'$\ell(\ell+1)P_{n\kappa}(\ell)/2\pi$')
	ax.set_title('W%i_zcut%shi_sigmaG%02d'%(Wx,zcut,sigmaG*10))
	#show()
	savefig(plot_dir+'CC_edges_W%i_zcut%shi_sigmaG%02d.jpg'%(Wx,zcut,sigmaG*10))
	close()

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def KSmap (i):
	'''
	Smooth and KS inversion
	input: i = 1, 2, 3 ... 13
	output: nothing, write kmap, Mmask if haven't done so
	'''
	Me1, Me2, Mw, galn = fileGen(i)
	for sigmaG in sigmaG_arr:	
		print 'KSmap i, sigmaG', i, sigmaG
		KS_fn = KS_dir+'CFHT_KS_sigma%02d_subfield%02d.fits'%(sigmaG*10,i)
		mask_fn = '/scratch/02977/jialiu/KSsim/mask/CFHT_mask_ngal%i_sigma%02d_subfield%02d.fits'%(ngal_arcmin,sigmaG*10,i)
		
		if WLanalysis.TestComplete((KS_fn,mask_fn),rm=True):
			kmap = WLanalysis.readFits(KS_fn)
			Mmask = WLanalysis.readFits(mask_fn)
		else:
			Me1_smooth = WLanalysis.weighted_smooth(Me1, Mw, PPA=PPA512, sigmaG=sigmaG)
			Me2_smooth = WLanalysis.weighted_smooth(Me2, Mw, PPA=PPA512, sigmaG=sigmaG)
			galn_smooth = snd.filters.gaussian_filter(galn.astype(float),sigmaG*PPA512, mode='constant')
			## KS
			kmap = WLanalysis.KSvw(Me1_smooth, Me2_smooth)
			## mask
			maskidx = where(galn_smooth < ngal_cut) #cut at ngal=5
			Mmask = ones(shape=galn.shape)
			Mmask[maskidx]=0
			
			WLanalysis.writeFits(kmap, KS_fn)
			WLanalysis.writeFits(Mmask, mask_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [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 readFits [as 别名]
	def ips_pk_single (R):#, sigmaG, zg, bins):
		kmap = WLanalysis.readFits(KSsim_fn(i, cosmo, R, sigmaG, zg))
		if pk:#peaks
			mask = WLanalysis.readFits(Mask_fn(i, sigmaG))
			peaks_hist = WLanalysis.peaks_mask_hist(kmap, mask, bins, kmin=kmin, kmax=kmax)
			return peaks_hist
		else:#powspec
			ell_arr, powspec = WLanalysis.PowerSpectrum(kmap, sizedeg=12.0)
			return powspec

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def average_powspec_nonoise (cosmo):
	ps = zeros(shape=(1000,50))
	weights = (genfromtxt(KSsim_dir+'galn.txt').T[1]).astype(float)
	weights /= sum(weights)
	fn = KSsim_dir+'powspec_Mk_sum13fields/SIM_powspec_sigma05_rz1_%s_1000R.fit'%(cosmo)
	if os.path.isfile(fn):
		return WLanalysis.readFits(fn)
	else:
		for i in range(1,14):
			ips=weights[i-1]*WLanalysis.readFits(powspecMk_fn(i, cosmo))
			ps += ips
		WLanalysis.writeFits(ps,fn)
		return ps

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
	def z_hist (i):
		'''return the histogramed redshift distribution for subfield i'''
		fn = '/Users/jia/weaklensing/CFHTLenS/catalogue/emulator_galpos_zcut0213/emulator_subfield%i_zcut0213.fit'%(i)
		z = WLanalysis.readFits(fn).T[2]
		zhist = histogram(z, range=(0.2,1.3), bins=16)
		print i, len(z)
		return zhist[0]

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
	def plotemupk (sigmabins):

		sigmaG, bins = sigmabins
		x = linspace(-0.04, 0.12, bins+1)
		x = x[:-1]+0.5*(x[1]-x[0])
		print "sigmaG, bins", sigmaG, bins
		def getpk (pk_fn, bins = bins):
			pk600bins = WLanalysis.readFits(emu_dir+'peaks_sum/sigma%02d/'%(sigmaG*10)+pk_fn)
			pk = pk600bins.reshape(1000, -1, 600/bins)
			pk = sum(pk, axis = -1)
			return pk
		CFHT_peak = zeros(bins)
		for j in arange(1,14):
			#print 'adding up subfield,',j
			CFHT_peak+=WLanalysis.readFits (CFHT_dir+'CFHT_peaks_sigma%02d_subfield%02d_%03dbins.fits'%(sigmaG*10, j, bins))
			
		pk_fn_arr = os.listdir(emu_dir+'peaks_sum/sigma%02d/'%(sigmaG*10))
		pk_mat = array(map(getpk, pk_fn_arr))
		pk_avg = mean(pk_mat,axis=1)
		pk_std = std(pk_mat, axis=1)
		
		for i in range(len(pk_avg)):
			plot(x, pk_avg[i], color=rand(3))
		plot(x, CFHT_peak, 'k--', linewidth = 2)
		xlabel('kappa')
		title('nbin %i, sigmaG %s'%(bins, sigmaG))
		#show()
		savefig(plot_dir+'emu_peaks_test_bins%i_sigmaG%02d.jpg'%(bins, sigmaG*10))
		close()

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def Psingle_CFHT (i, sigmaG, bins, ps=0):
	if ps:
		fn = powspec_CFHT_fn(i, sigmaG)
	else:
		fn = peaks_CFHT_fn(i, sigmaG, bins)
	
	if os.path.isfile(fn):
		out = WLanalysis.readFits(fn)
	elif ps:
		kmap = WLanalysis.readFits(KSCFHT_fn(i, sigmaG))
		out = WLanalysis.PowerSpectrum(kmap, sizedeg=12.0)
		WLanalysis.writeFits(out,fn)
	else:
		kmap = WLanalysis.readFits(KSCFHT_fn(i, sigmaG))
		mask = WLanalysis.readFits(Mask_fn(i, sigmaG))
		out = WLanalysis.peaks_mask_hist(kmap, mask, bins, kmin=kmin, kmax=kmax)
		WLanalysis.writeFits(out,fn)
	return out

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def fileGen(i):
	'''
	Input:
	i range from (1, 2..13)
	Return:
	Me1 = e1*w
	Me2 = (e2-c2)*w
	Mw = (1+m)*w
	galn = number of galaxies per pixel
	'''
	Me1_fn = KS_dir+'CFHT_subfield%02d_Me1.fits'%(i)
	Me2_fn = KS_dir+'CFHT_subfield%02d_Me2.fits'%(i)
	Mw_fn = KS_dir+'CFHT_subfield%02d_Mw.fits'%(i)
	galn_fn = KS_dir+'CFHT_subfield%02d_galn.fits'%(i)
	
	print 'fileGen', i
	if WLanalysis.TestComplete((Me1_fn,Me2_fn,Mw_fn,galn_fn),rm = True):
		Me1 = WLanalysis.readFits(Me1_fn)
		Me2 = WLanalysis.readFits(Me2_fn)
		Mw =  WLanalysis.readFits(Mw_fn)
		galn =WLanalysis.readFits(galn_fn)
	else:
		ifile = np.genfromtxt(full_dir+'full_subfield'+str(i) ,usecols=[0, 1, 2, 3, 4, 9, 10, 11, 16, 17])
		# cols: y, x, z_peak, z_rnd1, z_rnd2, e1, e2, w, m, c2

		#redshift cut 0.2< z <1.3
		zs = ifile[:,[2,3,4]]
		print 'zs'
		idx = np.where((amax(zs,axis=1) <= zmax) & (amin(zs,axis=1) >= zmin))[0]
		
		y, x, z_peak, z_rnd1, z_rnd2, e1, e2, w, m, c2 = ifile[idx].T

		k = array([e1*w, (e2-c2)*w, (1+m)*w])
		Ms, galn = WLanalysis.coords2grid(x, y, k)
		print 'coords2grid'
		Me1, Me2, Mw = Ms
		WLanalysis.writeFits(Me1,Me1_fn)
		WLanalysis.writeFits(Me2,Me2_fn)
		WLanalysis.writeFits(Mw,Mw_fn)
		WLanalysis.writeFits(galn,galn_fn)
	return Me1, Me2, Mw, galn

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def average_powspec_withnoise (cosmo, sigmaG, zg='rz1', CFHT=None):
	weights = (genfromtxt(KSsim_dir+'galn.txt').T[1]).astype(float)
	weights /= sum(weights)
	
	if CFHT:
		ps = zeros(50)
		fn = KSsim_dir+'powspec_sum13fields/CFHT_powspec_sigma%02d.fit'%(sigmaG*10)
	else:
		ps = zeros(shape=(1000,50))
		fn = KSsim_dir+'powspec_sum13fields/SIM_powspec_sigma%02d_%s_%s_%04dR.fit'%(sigmaG*10, zg, cosmo, Rtol)
	if os.path.isfile(fn):
		return WLanalysis.readFits(fn)
	else:
		for i in range(1,14):
			if CFHT:
				ips=weights[i-1]*WLanalysis.readFits(powspec_CFHT_fn(i, sigmaG))[-1]
			else:
				ips=weights[i-1]*WLanalysis.readFits(powspec_fn(i, cosmo, 1000, sigmaG, zg))
			ps += ips
		WLanalysis.writeFits(ps,fn)
		return ps

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def Noise(iinput):
	'''Input: (Wx, iseed)
	Return: files of noise KS map, using randomly rotated galaxy.
	'''
	Wx, iseed = iinput
	seed(iseed)
	print 'Bmode - Wx, iseed:', Wx, iseed
	bmap_fn = cat_dir+'Noise/W%i/W%i_Noise_sigmaG10_%04d.fit'%(Wx, Wx, iseed)
	
	isfile_kmap, bmap = WLanalysis.TestFitsComplete(bmap_fn, return_file = True)
	if isfile_kmap == False:
		Me1_fn = cat_dir+'Me_Mw_galn/W%i_Me1w_1.3_lo.fit'%(Wx)
		Me2_fn = cat_dir+'Me_Mw_galn/W%i_Me2w_1.3_lo.fit'%(Wx)
		Mw_fn = cat_dir+'Me_Mw_galn/W%i_Mwm_1.3_lo.fit'%(Wx)
		Me1_init = WLanalysis.readFits(Me1_fn)
		Me2_init = WLanalysis.readFits(Me2_fn)
		#### randomly rotate Me1, Me2 ###
		Me1, Me2 = WLanalysis.rndrot(Me1_init, Me2_init)
		#################################
		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)
		bmap = WLanalysis.KSvw(Me1_smooth, Me2_smooth)
		WLanalysis.writeFits(bmap,bmap_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
		def returnCC (Wx):
			edges = edgesGen(Wx)
			#print Wx
			mask = maskGen(Wx, 0.5, sigmaG)
			galn = WLanalysis.readFits(obsPK_dir+'maps/W%i_galn_1.3_lo_sigmaG%02d.fit'%(Wx, sigmaG*10))
			kmap = kmap_lensing_Gen(Wx, sigmaG)
			bmode = bmode_lensing_Gen(Wx, sigmaG)
			kproj = kmap_predict_Gen(Wx, sigmaG)
			
			ell_arr, pk = WLanalysis.CrossCorrelate(kmap*mask, galn*mask,edges=edges)
			ell_arr, pb = WLanalysis.CrossCorrelate(kmap*mask, bmode*mask,edges=edges)
			ell_arr, pp = WLanalysis.CrossCorrelate(kproj*mask, galn*mask,edges=edges)

			ell_arr, ppk = WLanalysis.CrossCorrelate(kproj*mask, kmap*mask,edges=edges)
			ell_arr, ppb = WLanalysis.CrossCorrelate(kproj*mask, bmode*mask,edges=edges)
			return ell_arr, pk, pb, pp, ppk, ppb

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def fileGen(i, R, cosmo):
	'''
	Put catalogue to grid, with (1+m)w correction. Mw is already done.
	also add randomly rotated noise
	Input:
	i: subfield range from (1, 2..13)
	R: realization range from (1..1000)
	cosmo: one of the 100 cosmos
	Return:
	Me1 = e1*w
	Me2 = e2*w
	
	'''
	#y, x, e1, e2, w, m = yxewm_arr[i-1].T
	s1, s2 = (WLanalysis.readFits(SIMfn(i,cosmo,R)).T)[[1,2]]
	A, galn = WLanalysis.coords2grid(x, y, array([s1*w, s2*w]))
	Ms1, Ms2 = A
	return Ms1, Ms2

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import readFits [as 别名]
def createBadFieldMask (sf):
	sf_splitfiles = os.listdir(sf_dir(sf))
	genfromtxtA = lambda fn: genfromtxt(sf_dir(sf)+fn)
	datas = map(genfromtxtA,sf_splitfiles)#3 columns: RA, DEC, GB
	datas = concatenate(datas,axis=0)
	idx = where(datas[:,-1]==1)[0]
	datas = datas[idx]
	y, x, k = datas.T
	k, galn = WLanalysis.coords2grid(x, y, array([k,]))
	for sigmaG in sigmaG_arr:
		print 'createBadFieldMask sf, sigmaG:', sf, sigmaG
		Allmask = WLanalysis.readFits(mask_fcn(sigmaG, sf))#mask for all field
		badmask_fn = badmask_fcn(sigmaG, sf)#file name for bad pointing mask, which is 75% area of Allmask
		galn_smooth = snd.filters.gaussian_filter(galn.astype(float),sigmaG*PPA512, mode='constant')
		#smooth the galn grid
		Mmask = ones(shape=galn.shape)#create mask grid
		Mmask[where(galn_smooth < ngal_cut)]=0#find the low density region in galn_smooth
		Mmask = adHocFix(Mmask,sf)
		Mmask *= Allmask#since I didn't do redshift cut in badmask, so here it takes care of it, since ALl mask has redshift cuts
		WLanalysis.writeFits(Mmask, badmask_fn)