Python WLanalysis.edge2center方法代码示例

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import edge2center [as 别名]
def FT_PowerSpectrum (cosmo, r, bins=10, return_ell_arr=0, Gaus=0):
    if Gaus:
        a = FTmapGen_Gaus(r)
    else:
        a = FTmapGen(cosmo, r)
    PS2D=np.abs(fftpack.fftshift(a))**2
    ell_arr,psd1D=WLanalysis.azimuthalAverage(PS2D, bins=bins)    
    if return_ell_arr:
        ell_arr = WLanalysis.edge2center(ell_arr)* 360./3.5
        return ell_arr
    else:
        return psd1D

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import edge2center [as 别名]
	center = np.array([(x.max()-x.min())/2.0, (x.max()-x.min())/2.0])
	if n%2 == 0:
		center+=0.5
	r = np.hypot(x - center[0], y - center[1])#distance to center pixel, for each pixel
	r_sorted = r.flat # the index to sort by r

	# find index that's corresponding to the lower edge of each bin
	kmin=1.0
	kmax=n/2.0
	edges = edgesGen(Wx)		
	hist_ind = np.histogram(r_sorted,bins = edges)[0]
	N_indep = hist_ind[1:]
	return N_indep/2

if compute_theory_err:
	ell_arr = 40.0*WLanalysis.edge2center(linspace(1,50,6))
	b_ell = exp(-ell_arr**2*radians(1.0/60)**2/2.0)
	factor = 2.0*pi/(ell_arr+1)
	d_ell = ell_arr[1]-ell_arr[0]
	edgesGen = lambda Wx: linspace(1,50,6)*sizes[Wx-1]/1330.0
	sizes = (1330, 800, 1120, 950)
	sizedeg_arr = array([(sizes[Wx-1]/512.0)**2*12.0 for Wx in range(1,5)])
	kmapGen = lambda Wx: np.load(cmb_dir+'cfht/kmap_W%i_sigma10_noZcut.npy'%(Wx))
	maskGen = lambda Wx: np.load(cmb_dir+'mask/W%i_mask1315_noZcut.npy'%(Wx))
	cmblGen = lambda Wx: np.load(cmb_dir+'planck/COM_CompMap_Lensing_2048_R1.10_kappa_CFHTLS_W%i.npy'%(Wx))#2013
	#cmblGen = lambda Wx: np.load(cmb_dir+'planck/dat_kmap_flipper2048_CFHTLS_W%i_map.npy'%(Wx))#2015
	edges_arr = map(edgesGen, range(1,5))
	mask_arr = map(maskGen, range(1,5))
	fmask_arr = array([sum(mask_arr[Wx-1])/sizes[Wx-1]**2 for Wx in range(1,5)])
	fmask2_arr = array([sum(mask_arr[Wx-1]**2)/sizes[Wx-1]**2 for Wx in range(1,5)])
	fsky_arr = fmask_arr*sizedeg_arr/41253.0

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import edge2center [as 别名]
	ax.plot(S_arr, Pps_marg,'g--', label=r'$\rm{power\,spectrum}\,(\alpha=0.64)$',linewidth=2)
	ax.plot(S_arr, Ppk_marg,'m-',label=r'$\rm{peaks}\,(\alpha=0.60)$',linewidth=1)
	ax.plot(S_arr, Pcomb_marg, 'k-',label=r'$\rm{power\, spectrum + peaks}\,(\alpha=0.63)$',linewidth=2)
	leg=ax.legend(ncol=1, labelspacing=0.3, prop={'size':16},loc=0)
	leg.get_frame().set_visible(False)
	ax.set_xlabel(r'$\rm{\Sigma_8=\sigma_8(\Omega_m/0.27)^\alpha}$',fontsize=20)
	ax.set_ylabel(r'$\rm{Probability}$',fontsize=20)
	ax.set_xlim(0.5, 1.2)
	ax.set_ylim(0.0, 0.07)
	#ax.set_title('w=%s'%(w_arr[i]))
	#savefig(plot_dir+'SIGMA_marg_prob_w%s.pdf'%(w_arr[i]))
	savefig(plot_dir+'SIGMA_marg_prob.pdf')
	close()

if sample_points:
	kappa_center=WLanalysis.edge2center(linspace(-.04,.12,26))
	ell_arr_labels = array(['%.3f'%(i) for i in list(kappa_center.copy())*2])
	best_fit_arr = ([0.44, -0.78, 0.62],[0.42, -0.87, 0.67],[0.30, -0.78, 0.82])
	
	import test_chisq_cube_MPI as tcs
	interp_cosmo, cov_mat, cov_inv, ps_CFHT = tcs.return_interp_cosmo_for_idx (tcs.idx_pk2)

	f=figure(figsize=(10,8))
	ax=f.add_subplot(gs[0])
	ax2=f.add_subplot(gs[1],sharex=ax)
	
	fidu_std = sqrt(diag(cov_mat))
	fidu_avg = tcs.ps_avg0[tcs.idx_pk2]
	
	ax.errorbar(ell_arr, ps_CFHT, fidu_std, color='k', linewidth=lw)
	ax.plot(ell_arr, ps_CFHT, 'k-', label=r'$\rm{CFHTLenS}$',linewidth=lw)	

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import edge2center [as 别名]
#CmaskGen = lambda Wx: load (main_dir+'cfht_mask/Mask_W%s_0.7_sigmaG10.npy'%(Wx))
CmaskGen = lambda Wx: load (main_dir+'mask_ludo/ludomask_weight0_manu_W%i.npy'%(Wx))
PmaskGen = lambda Wx: load (main_dir+'planck2015_mask/kappamask_flipper2048_CFHTLS_W%s_map.npy'%(Wx))
maskGen = lambda Wx: CmaskGen(Wx)*PmaskGen(Wx)
PlanckSim15Gen = lambda Wx, r: load('/work/02977/jialiu/cmblensing/planck/sim15/sim_%04d_kmap_CFHTLS_W%s.npy'%(r, Wx))



edgesGen = lambda Wx: linspace(1,60,7)*sizes[Wx-1]/1330.0
### omori & holder bin edges #####
#edgesGen = lambda Wx: linspace(1.25, 47.49232195,21)*sizes[Wx-1]/1330.0

edges_arr = map(edgesGen, range(1,5))
sizedeg_arr = array([(sizes[Wx-1]/512.0)**2*12.0 for Wx in range(1,5)])
####### test: ell_arr = WLanalysis.PowerSpectrum(CmaskGen(1), sizedeg = sizedeg_arr[0],edges=edges_arr[0])[0]
ell_arr = WLanalysis.edge2center(edgesGen(1))*360.0/sqrt(sizedeg_arr[0])

factor = (ell_arr+1)*ell_arr/(2.0*pi)
mask_arr = map(maskGen, range(1,5))
fmask_arr = array([sum(mask_arr[Wx-1])/sizes[Wx-1]**2 for Wx in range(1,5)])
fmask2_arr = array([sum(mask_arr[Wx-1]**2)/sizes[Wx-1]**2 for Wx in range(1,5)])
fsky_arr = fmask_arr*sizedeg_arr/41253.0
d_ell = ell_arr[1]-ell_arr[0]
#################################################

def theory_CC_err(map1, map2, Wx):
    map1*=mask_arr[Wx-1]
    map2*=mask_arr[Wx-1]
    #map1-=mean(map1)
    #map2-=mean(map2)
    auto1 = WLanalysis.PowerSpectrum(map1, sizedeg = sizedeg_arr[Wx-1], edges=edges_arr[Wx-1],sigmaG=1.0)[-1]/fmask2_arr[Wx-1]/factor

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import edge2center [as 别名]
    for iom, isi8 in cosmo_params_noisy:
        if not iom==om0:
            ax.scatter(iom,isi8, s=300, marker='o',color='deepskyblue', facecolors='none', edgecolors='deepskyblue',lw=2)
    ax.set_xlabel(r'$\Omega_m$',fontsize=22)
    ax.set_ylabel(r'$\sigma_8$',fontsize=22)
    ax.set_xlim(0.12,0.75)
    ax.set_ylim(0.46,1.05)
    ax.grid(True)
    #show()
    savefig(CMBNG_dir+'plot_official/plot_design.pdf')
    savefig(CMBNG_dir+'plot/plot_design.png')
    close()
    
if plot_comp_nicaea:
    cosmo='Om0.296_Ol0.704_w-1.000_si0.786'
    ell_gadget = (WLanalysis.edge2center(logspace(log10(1.0),log10(1024),51))*360./sqrt(12.25))[:34]

    ell_nicaea, ps_nicaea=genfromtxt('/Users/jia/weaklensing/CMBnonGaussian/Pkappa_nicaea/Pkappa_nicaea25_{0}_1100'.format(cosmo))[33:-5].T
    
    ell_nicaea2, ps_nicaea_linear=genfromtxt('/Users/jia/weaklensing/CMBnonGaussian/Pkappa_nicaea/Pkappa_nicaea25_{0}_1100_linear'.format(cosmo))[33:-5].T
    
    def get_1024(j):
        pspkPDFgadget=load('/Users/jia/weaklensing/CMBnonGaussian/Pkappa_gadget/noiseless/kappa_{0}_ps_PDF_pk_z1100_{1}.npy'.format(fidu_cosmo, j))
        ps_gadget=array([pspkPDFgadget[i][0] for i in range(len(pspkPDFgadget))])
        ps_gadget=ps_gadget.squeeze()
        return ps_gadget
    ps_gadget=concatenate(map(get_1024, range(10)),axis=0)[:,:34]
    
    idx=where(~isnan(mean(ps_gadget,axis=0)))[0]
    ps_gadget=ps_gadget[:,idx]
    ell_gadget=ell_gadget[idx]

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import edge2center [as 别名]
	all_kappa = kmaps[~isnan(kmaps)]
	PDF = histogram(all_kappa, bins=edges)[0]
	PDF_normed = PDF/float(len(all_kappa))
	return PDF_normed

#pool = MPIPool()
#for cosmo in cosmo_arr:
	#print cosmo
	#if not os.path.isfile(fn(cosmo)):
		#cosmoR_arr = [(cosmo, R) for R in range(1,1001)]
		#PDF_arr = array(pool.map(PDFGen, cosmoR_arr))
		#np.save(fn(cosmo), PDF_arr)

#fsky_all = array([0.839298248291,0.865875244141,0.809467315674,
		  #0.864688873291,0.679264068604,0.756385803223,
		  #0.765892028809,0.747268676758,0.77250289917,
		  #0.761451721191,0.691867828369,0.711254119873,
		  #0.745429992676])

PDF92 = mean(array([load(fn(cosmo)) for cosmo in cosmo_arr]), axis=1)
for iPDF in PDF92:
	plot(WLanalysis.edge2center(edges), iPDF)
if noise:
	title('noisy')
else:
	title('noiseless')
xlabel('kappa')
ylabel('PDF')
	
show()
print 'done done done'