Python WLanalysis.weighted_smooth方法代码示例

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [as 别名]
def KSmap_massproduce(iiRcosmo):
	'''Input:
	iiRcosmo = [i, R, cosmo]
	i: subfield range from (1, 2..13)
	R: realization range from (1..1000)
	cosmo: one of the 1000 cosmos
	Return:
	KS inverted map
	'''
	i, R, cosmo = iiRcosmo
	create_kmap = 0
	for sigmaG in sigmaG_arr:
		KS_fn = KSfn(i, cosmo, R, sigmaG)
		if not os.path.isfile(KS_fn):
			create_kmap = 1
			break
	if create_kmap:
		print 'Mass Produce, creating KS: ', i, R, cosmo
		Me1, Me2 = fileGen(i, R, cosmo)
		for sigmaG in sigmaG_arr:
			KS_fn = KSfn(i, cosmo, R, sigmaG)	
			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)
			np.save(KS_fn, kmap)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_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 weighted_smooth [as 别名]
	def randmap (iseed, Wx=Wx):	
		Me1rnd, Me2rnd = WLanalysis.rndrot(Me1, Me2, iseed=iseed)
		Me1smooth = WLanalysis.weighted_smooth(Me1rnd, Mwm)
		Me2smooth = WLanalysis.weighted_smooth(Me2rnd, Mwm)
		kmap_rand = WLanalysis.KSvw(Me1smooth, Me2smooth)
		print Wx, iseed, kmap_rand.shape
		np.save(bmap_fn(Wx, iseed), kmap_rand)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [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 weighted_smooth [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 weighted_smooth [as 别名]
	def randmap (iseedWx):
		iseed, Wx = iseedWx
		Me1, Me2 = Me1_arr[Wx-1], Me2_arr[Wx-1]
		Mwm = Mwm_arr[Wx-1]
		Me1rnd, Me2rnd = WLanalysis.rndrot(Me1, Me2, iseed=iseed)
		Me1smooth = WLanalysis.weighted_smooth(Me1rnd, Mwm)
		Me2smooth = WLanalysis.weighted_smooth(Me2rnd, Mwm)
		kmap_rand = WLanalysis.KSvw(Me1smooth, Me2smooth)
		print Wx, iseed, kmap_rand.shape
		np.save(bmap_fn(Wx, iseed), kmap_rand)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [as 别名]
def kmapPk_1sim (r):
	print r,'1sim'
	k, s1, s2 = kappaGen_1sim(r)[:3]
	A, galn = WLanalysis.coords2grid(x, y, array([s1,s2]))
	Ms1,Ms2 = A
	s1_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	s2_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	kmap = WLanalysis.KSvw(s1_smooth, s2_smooth)
	pk = WLanalysis.peaks_mask_hist(kmap, mask, bins=25, kmin = -0.04, kmax = 0.12)
	return pk

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [as 别名]
def kmapPs (r):
	print r
	k, s1, s2 = kappaGen(r)[:3]
	A, galn = WLanalysis.coords2grid(x, y, array([s1,s2 ]))
	Ms1,Ms2 = A
	s1_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	s2_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	kmap = WLanalysis.KSvw(s1_smooth, s2_smooth)
	kmap*=mask
	ps = WLanalysis.PowerSpectrum(kmap,sizedeg=12.0)[-1]
	return ps

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [as 别名]
def KSmap_single(i, R, cosmo, sigmaG):
	'''Input:
	i: subfield range from (1, 2..13)
	R: realization range from (1..1000)
	cosmo: one of the 1000 cosmos
	sigmaG: smoothing scale
	Return:
	KS inverted map
	'''
	KS_fn = KSfn(i, cosmo, R, sigmaG)
	if not os.path.isfile(KS_fn):
		Me1, Me2 = fileGen(i, R, cosmo)
		#Mw = Mw_fcn(i)
		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)
		np.save(KS_fn, kmap)
	return kmap

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [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 weighted_smooth [as 别名]
		for sG in sG_arr:
			KS_45_fn = test_dir+'KS_45_%i_%i.fit'%(i,sG)
			KS_rand_fn = test_dir+'KS_rand_%i_%i.fit'%(i,sG)
			print i, sG
			if os.path.isfile(KS_45_fn) and os.path.isfile(KS_rand_fn):
				KS_45=WLanalysis.readFits(KS_45_fn)
				KS_rand=WLanalysis.readFits(KS_rand_fn)
			
			else:
				print 'generating KS'
				y, x, e1, e2, w = WLanalysis.readFits(test_dir+'yxew_subfield%i_zcut0213.fit'%(i)).T
				e1_45, e2_45 = WLanalysis.rndrot(e1, e2, deg=45)
				e1_rand, e2_rand = WLanalysis.rndrot(e1, e2, iseed=0)
				mat_e1_45,mat_e2_45,mat_e1_rand,mat_e2_rand, Mw = WLanalysis.coords2grid(x, y, array([e1_45*w, e2_45*w, e1_rand*w, e2_rand*w, w]) )[0]

				mat_e1_45_smoothed  = WLanalysis.weighted_smooth(mat_e1_45  , Mw, sigmaG=sG)
				mat_e2_45_smoothed  = WLanalysis.weighted_smooth(mat_e2_45  , Mw, sigmaG=sG)
				mat_e1_rand_smoothed= WLanalysis.weighted_smooth(mat_e1_rand, Mw, sigmaG=sG)
				mat_e2_rand_smoothed= WLanalysis.weighted_smooth(mat_e2_rand, Mw, sigmaG=sG)
				KS_45 =WLanalysis.KSvw(mat_e1_45_smoothed,mat_e2_45_smoothed)
				KS_rand=WLanalysis.KSvw(mat_e1_rand_smoothed,mat_e2_rand_smoothed)
				
				
				WLanalysis.writeFits(KS_45,KS_45_fn)
				WLanalysis.writeFits(KS_rand,KS_rand_fn)
			
			kappa_45 = KS_45.flatten()
			std_45=std(kappa_45)
			#hist_45,binedges45 = histogram(kappa_45,bins=100,range=(-5*std_45,5*std_45))
			#savetxt(test_dir+'hist45_%i_%i.ls'%(i,sG),array([hist_45,binedges45[:-1]]).T)
			

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [as 别名]
		#####################################
		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)
		
		mask_nan = mask0.copy()
		mask_nan[mask0==0]=nan

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [as 别名]
def ps (R):
	Mk = WLanalysis.readFits(Mk_fn(i, cosmo, R))
	galn = WLanalysis.readFits(galn_fn(i))
	Mk_smooth = WLanalysis.weighted_smooth(Mk, galn, sigmaG)
	ps = WLanalysis.PowerSpectrum(Mk_smooth)
	return ps[1]

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import weighted_smooth [as 别名]
	s2 = s2o*(1+m)
	
	eint1, eint2 = WLanalysis.rndrot(e1, e2, iseed=1000)
	eint1_45, eint2_45 = WLanalysis.rndrot(e1, e2, iseed=1000,deg=45.0)

	e1_reduce, e2_reduce = eobs_fun(s1, s2, kappa, eint1, eint2)
	e1_add, e2_add = s1+eint1, s2+eint2
	
	A, galn = WLanalysis.coords2grid(x, y, array([e1_reduce*w, e2_reduce*w, w*(1+m), s1o, s2o, kappa, e1_add*w, e2_add*w, s1+znoise1, s2+znoise2, kappa+znoise1]))
	#Mk, Ms1, Ms2 = A
	Me1, Me2, Mw, Ms1, Ms2, Mk, Me1add, Me2add, Ms1n, Ms2n, Mkn = A
	
	B, galn = WLanalysis.coords2grid(x, y, array([znoise1, eint1, eint2, eint1_45, eint2_45, znoise27]))
	Mn, Men1, Men2, Men1b, Men2b, Mn27 = B
	### pure analytical noise ######
	Mn_smooth =  WLanalysis.weighted_smooth(Mn, galn, sigmaG=0.5)
	Mn27_smooth =  WLanalysis.weighted_smooth(Mn27, galn, sigmaG=0.5)
	Men1_smooth =  WLanalysis.weighted_smooth(Men1, galn, sigmaG=0.5)
	Men2_smooth =  WLanalysis.weighted_smooth(Men2, galn, sigmaG=0.5)
	Men1b_smooth =  WLanalysis.weighted_smooth(Men1b, galn, sigmaG=0.5)
	Men2b_smooth =  WLanalysis.weighted_smooth(Men2b, galn, sigmaG=0.5)
	Nmap = WLanalysis.KSvw(Men1_smooth, Men2_smooth)
	Nmap45 = WLanalysis.KSvw(Men1b_smooth, Men2b_smooth)
	
	ps_nconv27 = WLanalysis.PowerSpectrum(Mn27_smooth)[-1]
	ps_nconv = WLanalysis.PowerSpectrum(Mn_smooth)[-1]
	ps_nrand = WLanalysis.PowerSpectrum(Nmap)[-1]
	ps_n45 = WLanalysis.PowerSpectrum(Nmap45)[-1]
	
	loglog(ell_arr, ps_nconv27, '.',label='Conv pure noise rms=0.29')
	loglog(ell_arr, ps_nconv, '.', label='Conv pure noise rms=0.15+0.035z')