From f12c9a0c1af0f997858f591238b11ca43f138e70 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 22 Jan 2015 19:16:46 +0100 Subject: [PATCH] KSZ work --- .../build_dipole_ksz_from_galaxies.py | 76 ++++++++++++++++--- python_sample/icgen/borgicgen.py | 10 ++- python_sample/icgen/cosmogrowth.py | 36 ++++++--- python_sample/ksz/constants.py | 1 + python_sample/ksz/gal_prof.py | 38 ++++++++-- 5 files changed, 128 insertions(+), 33 deletions(-) diff --git a/python_sample/build_dipole_ksz_from_galaxies.py b/python_sample/build_dipole_ksz_from_galaxies.py index 77eca95..2d1d200 100644 --- a/python_sample/build_dipole_ksz_from_galaxies.py +++ b/python_sample/build_dipole_ksz_from_galaxies.py @@ -100,13 +100,35 @@ def build_unit_vectors(N): return ux,uy,uz -def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_hubble=False): +def compute_vcmb(l, b): + # Motion is obtained from Tully (2007): sun_vs_LS + LS_vs_CMB + motion = [-25.,-246.,277.]; + + x = np.cos(l*np.pi/180) * np.cos(b*np.pi/180) + y = np.sin(l*np.pi/180) * np.cos(b*np.pi/180) + z = np.sin(b*np.pi/180) + + return x*motion[0] + y*motion[1] + z*motion[2] + + +def compute_vlg(l,b): + + motion = [-79,296,-36]; # [-86, 305, -33]; + + x = np.cos(l*np.pi/180) * np.cos(b*np.pi/180) + y = np.sin(l*np.pi/180) * np.cos(b*np.pi/180) + z = np.sin(b*np.pi/180) + + return x*motion[0] + y*motion[1] + z*motion[2] + + +def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_hubble=False,x=2.37,bright=-np.inf,bright_list=[],use_vlg=True,sculpt=-1): import progressbar as pbar cat = np.load("2m++.npy") cat['distance'] = cat['best_velcmb'] - cat = cat[np.where((cat['distance']>100*dmin)*(cat['distance']100*dmin)*(cat['distance'] dmax: + continue + Lgal = DA**2*10**(0.4*(tmpp_cat['Msun']-i['K2MRS']+25)) - profiler = ksz.KSZ_Isothermal(Lgal, 2.37, y=y) + M_K=i['K2MRS']-5*np.log10(DA)-25 + # Skip too bright galaxies + if M_K < bright: + continue + + profiler = ksz.KSZ_Isothermal(Lgal, x, y=y, sculpt=sculpt) idx0 = hp.query_disc(Nside, (x0,y0,z0), 3*profiler.rGalaxy/DA) @@ -177,12 +217,17 @@ def get_args(): parser.add_argument('--depth_max', type=float, default=60) parser.add_argument('--ksz_map', type=str, required=True) parser.add_argument('--base_fig', type=str, default="kszfig.png") - parser.add_argument('--build_dipole', type=bool, default=False) + parser.add_argument('--build_dipole', action='store_true') parser.add_argument('--degrade', type=int, default=-1) parser.add_argument('--y',type=float,default=0.0) - parser.add_argument('--random', type=bool, default=False) + parser.add_argument('--x',type=float,default=2.37) + parser.add_argument('--random', action='store_true') parser.add_argument('--perturb', type=float, default=0) - parser.add_argument('--hubble_monopole', type=bool, default=False) + parser.add_argument('--hubble_monopole', action='store_true') + parser.add_argument('--remove_bright', type=float, default=-np.inf) + parser.add_argument('--bright_file', type=str) + parser.add_argument('--lg', action='store_true') + parser.add_argument('--sculpt_beam', type=float, default=-1) return parser.parse_args() def main(): @@ -195,7 +240,16 @@ def main(): print("Generating map...") - r = generate_from_catalog(args.depth_min,args.depth_max,args.Nside,perturb=args.perturb,y=args.y,do_random=args.random,do_hubble=args.hubble_monopole) + with open("crap.txt", mode="r") as f: + bright_list = [l.split('#')[0].strip(" \t\n\r") for l in f] + + if args.bright_file: + with open(args.bright_file, mode="r") as f: + idx_name = f.readline().split(',').index('name_2') + bright_list = bright_list + [l.split(',')[idx_name] for l in f] + + print("Built bright point source list: " + repr(bright_list)) + r = generate_from_catalog(args.depth_min,args.depth_max,args.Nside,perturb=args.perturb,y=args.y,do_random=args.random,do_hubble=args.hubble_monopole,x=args.x,bright=args.remove_bright,use_vlg=args.lg,bright_list=bright_list,sculpt=args.sculpt_beam) hubble_map = None if args.hubble_monopole: proj,mask,hubble_map = r diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index fcfee2c..f5653b3 100644 --- a/python_sample/icgen/borgicgen.py +++ b/python_sample/icgen/borgicgen.py @@ -4,10 +4,11 @@ import cosmolopy as cpy from cosmogrowth import * import borgadaptor as ba -def gen_posgrid(N, L): +@ct.timeit +def gen_posgrid(N, L, delta=1, dtype=np.float32): """ Generate an ordered lagrangian grid""" - ix = (np.arange(N)*L/N).astype(np.float32) + ix = (np.arange(N)*(L/N*delta)).astype(dtype) x = ix[:,None,None].repeat(N, axis=1).repeat(N, axis=2) y = ix[None,:,None].repeat(N, axis=0).repeat(N, axis=2) @@ -147,6 +148,8 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, supergenerate D2 = -3./7 * D1_0**2 + if do_lpt2: + psi2 = lpt.lpt2('all') for j in xrange(3): # Generate psi_j (displacement along j) print("LPT1 axis=%d" % j) @@ -154,8 +157,7 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, supergenerate psi = psi.reshape((psi.size,)) if do_lpt2: print("LPT2 axis=%d" % j) - psi2 = lpt.lpt2(j) - psi += D2 * psi2.reshape((psi2.size,)) + psi += D2 * psi2[j].reshape((psi2[j].size,)) # Generate posx posx.append(((posq[j] + psi)%L).astype(np.float32)) # Generate vel diff --git a/python_sample/icgen/cosmogrowth.py b/python_sample/icgen/cosmogrowth.py index 234524c..f5f3652 100644 --- a/python_sample/icgen/cosmogrowth.py +++ b/python_sample/icgen/cosmogrowth.py @@ -4,6 +4,7 @@ import pyfftw import weakref import numpy as np import cosmolopy as cpy +import cosmotool as ct class CubeFT(object): def __init__(self, L, N, max_cpu=-1): @@ -98,6 +99,7 @@ class LagrangianPerturbation(object): self._kz = self.ik[None,None,:(self.N/2+1)] self.cache = {}#weakref.WeakValueDictionary() + @ct.timeit_quiet def upgrade_sampling(self, supersample): N2 = self.N * supersample N = self.N @@ -113,14 +115,26 @@ class LagrangianPerturbation(object): self.N = N2 self.cube = CubeFT(self.L, self.N, max_cpu=self.max_cpu) + @ct.timeit_quiet def _gradient(self, phi, direction): - ne.evaluate('phi_hat * i * kv / (kx**2 + ky**2 + kz**2)', out=self.cube.dhat, - local_dict={'i':-1j, 'phi_hat':phi, 'kv':self._kdir(direction), + if direction == 'all': + dirs = [0,1,2] + copy = True + else: + dirs = [direction] + copy = False + ret=[] + for dir in dirs: + ne.evaluate('phi_hat * i * kv / (kx**2 + ky**2 + kz**2)', out=self.cube.dhat, + local_dict={'i':-1j, 'phi_hat':phi, 'kv':self._kdir(dir), 'kx':self._kx, 'ky':self._ky, 'kz':self._kz},casting='unsafe') # self.cube.dhat = self._kdir(direction)*1j*phi - self.cube.dhat[0,0,0] = 0 - return self.cube.irfft() + self.cube.dhat[0,0,0] = 0 + x = self.cube.irfft() + ret.append(x.copy() if copy else x) + return ret[0] if len(ret)==1 else ret + @ct.timeit_quiet def lpt1(self, direction=0): return self._gradient(self.dhat, direction) @@ -155,6 +169,7 @@ class LagrangianPerturbation(object): self.cube.density = array return self.cube.rfft() + @ct.timeit_quiet def lpt2(self, direction=0): # k2 = self._get_k2() # k2[0,0,0] = 1 @@ -170,12 +185,13 @@ class LagrangianPerturbation(object): for j in xrange(3): q = self._do_irfft( potgen0(j) ).copy() for i in xrange(j+1, 3): - ne.evaluate('div + q * pot', out=div_phi2, - local_dict={'div':div_phi2, 'q':q,'pot':self._do_irfft( potgen0(i), copy=False ) } - ) - ne.evaluate('div - pot**2',out=div_phi2, - local_dict={'div':div_phi2,'pot':self._do_irfft(potgen(i,j), copy=False) } - ) + with ct.time_block("LPT2 elemental (%d,%d)" %(i,j)): + ne.evaluate('div + q * pot', out=div_phi2, + local_dict={'div':div_phi2, 'q':q,'pot':self._do_irfft( potgen0(i), copy=False ) } + ) + ne.evaluate('div - pot**2',out=div_phi2, + local_dict={'div':div_phi2,'pot':self._do_irfft(potgen(i,j), copy=False) } + ) phi2_hat = self._do_rfft(div_phi2) #self.cache['lpt2_potential'] = phi2_hat diff --git a/python_sample/ksz/constants.py b/python_sample/ksz/constants.py index 558ad4d..2a76b27 100644 --- a/python_sample/ksz/constants.py +++ b/python_sample/ksz/constants.py @@ -27,6 +27,7 @@ tmpp_cat={'Msun':3.29, baryon_fraction = Omega_baryon / Omega_matter ksz_normalization = -T_cmb*sigmaT*v_unit/(lightspeed*mu*mp) * baryon_fraction +assert ksz_normalization < 0 rho_mean_matter = Omega_matter * (3*(100e3/Mpc)**2/(8*np.pi*G)) Lbar = tmpp_cat['lbar'] / Mpc**3 M_over_L_galaxy = rho_mean_matter / Lbar diff --git a/python_sample/ksz/gal_prof.py b/python_sample/ksz/gal_prof.py index de52eba..b32e658 100644 --- a/python_sample/ksz/gal_prof.py +++ b/python_sample/ksz/gal_prof.py @@ -10,7 +10,14 @@ class KSZ_Profile(object): R_star= 0.0 # 15 kpc/h L_gal0 = 10**(0.4*(tmpp_cat['Msun']-tmpp_cat['Mstar'])) - def __init__(self): + def __init__(self,sculpt): + """Base class for KSZ profiles + + Arguments: + sculpt (float): If negative, do not sculpt. If positive, there will be a 2d + suppression of the profile with a radius given by sculpt (in arcmins). + """ + self.sculpt = sculpt * np.pi/180/60. self.rGalaxy = 1.0 def evaluate_profile(self, r): @@ -37,6 +44,11 @@ class KSZ_Profile(object): if tan_theta_2.size > 0: idx_mask = np.append(idx_mask,idx[tan_theta_2.argmin()]) + if self.sculpt > 0: + theta = np.arctan(tan_theta) + cond = theta < self.sculpt + m[cond] *= (theta[cond]/self.sculpt)**2 + return idx,idx_mask,m @@ -48,10 +60,20 @@ class KSZ_Isothermal(KSZ_Profile): sigma_FP=160e3 #m/s R_innergal = 0.030 - def __init__(self, Lgal, x, y=0.0): - "Support for Isothermal profile" + def __init__(self, Lgal, x, y=0.0, sculpt=-1): + """Support for Isothermal profile + + Arguments: + Lgal (float): Galaxy luminosity in solar units + x (float): extent of halo in virial radius units + + Keyword arguments: + y (float): Inner part where there is no halo + sculpt (float): If negative, do not sculpt. If positive, there will be a 2d + suppression of the profile with a radius given by sculpt (in arcmins). + """ - super(KSZ_Isothermal,self).__init__() + super(KSZ_Isothermal,self).__init__(sculpt) self.R_gal = 0.226 * x self.R_innergal *= y @@ -71,20 +93,20 @@ class KSZ_Isothermal(KSZ_Profile): Q = np.zeros(r.size) - cond = (r<=0) + cond = (r<=1e-10) Q[cond] = rho0*2/Mpc * (rGalaxy-rInner)/(rGalaxy*rInner) cond = (r>0)*(r <= rInner) D['r'] = r[cond] Q[cond] = ne.evaluate('rho0*2/(Mpc*r) * (arctan(sqrt( (rGalaxy/r)**2 -1 )) - arctan(sqrt( (rInner/r)**2 - 1 )))', local_dict=D) - + cond = (r > rInner)*(r <= rGalaxy) D['r'] = r[cond] Q[cond] = ne.evaluate('rho0*2/(Mpc*r) * arctan(sqrt( (rGalaxy/r)**2 -1 ))', local_dict=D) - - return Q,np.where(r