Merge branch 'master' of bitbucket.org:glavaux/cosmotool

CMakeLists.txt

@@ -22,11 +22,11 @@ option(INTERNAL_HDF5 "Build internal version of HDF5" OFF)
 option(INTERNAL_NETCDF "Build internal version of NETCDF" OFF)
 option(INTERNAL_BOOST "Build internal version of BOOST" OFF)
 
-IF(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS)
+#IF(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS)
-  SET(CosmoTool_local CosmoTool_static)
+#  SET(CosmoTool_local CosmoTool_static)
-ELSE(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS)
+#ELSE(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS)
   SET(CosmoTool_local CosmoTool)
-ENDIF(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS)
+#ENDIF(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS)
 
 MESSAGE(STATUS "Using the target ${CosmoTool_local} to build python module")
 

python/cosmotool/simu.py

@@ -23,16 +23,16 @@ class SimulationBare(PySimulationBase):
   def merge(self, other):
 
     def _safe_merge(a, b):
-      if b:
+      if b is not None:
-        if a:
+        if a is not None:
           a = [np.append(q, r) for q,r in zip(a,b)]
         else:
           a = b
       return a
       
     def _safe_merge0(a, b):
-      if b:
+      if b is not None:
-        if a:
+        if a is not None:
           a = np.append(a, b)
         else:
           a = 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']<dmax*100))]
+#  cat = cat[np.where((cat['distance']>100*dmin)*(cat['distance']<dmax*100))]
 
   deg2rad = np.pi/180
   Npix = 12*Nside**2
@@ -119,14 +141,19 @@ def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_h
     ksz_hubble_template = np.zeros(ksz_template.size, dtype=np.float64)
 
   for i in pbar.ProgressBar(maxval = cat.size, widgets=[pbar.Bar(), pbar.ETA()])(cat):
+    # Skip too point sources
+    if i['name'] in bright_list:
+      print("Object %s is in bright list" % i['name'])
+      continue
+
     if do_random:
       l = np.random.rand()*360
       b = np.arcsin(2*np.random.rand()-1)*180/np.pi
     else:
-      l,b=i['gal_long'],i['gal_lat']
+      l0,b0=i['gal_long'],i['gal_lat']
 
-    l=ne.evaluate('l*deg2rad')
+    l=ne.evaluate('l0*deg2rad')
-    b=ne.evaluate('b*deg2rad')
+    b=ne.evaluate('b0*deg2rad')
 
     dtheta,dphi = np.random.randn(2)*perturb
     theta,l=move_direction(dtheta,dphi,0.5*np.pi - b, l)
@@ -137,10 +164,23 @@ def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_h
     y0 = np.sin(l)*np.cos(b)
     z0 = np.sin(b)
 
-    DA =i['distance']/100
+    if use_vlg:
+      vlg = i['best_velcmb'] - compute_vcmb(l0, b0) + compute_vlg(l0, b0)
+      DA = vlg/100
+    else:
+      DA = i['best_velcmb'] / 100
+
+    if DA < dmin or DA > 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

python_sample/icgen/borgicgen.py

@@ -4,10 +4,11 @@ import cosmolopy as cpy
 from cosmogrowth import *
 import borgadaptor as ba
 
+@ct.timeit
 def gen_posgrid(N, L, delta=1, dtype=np.float32):
     """ Generate an ordered lagrangian grid"""
 
-    ix = (np.arange(N)*L/N*delta).astype(dtype)
+    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[j].reshape((psi2[j].size,))
-          psi += D2 * psi2.reshape((psi2.size,))
         # Generate posx
         posx.append(((posq[j] + psi)%L).astype(np.float32))
         # Generate vel

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):
+        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(direction),
+              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()
+          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,6 +185,7 @@ class LagrangianPerturbation(object):
             for j in xrange(3):
                 q = self._do_irfft( potgen0(j) ).copy()
                 for i in xrange(j+1, 3):
+                    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 ) } 
                             )

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

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):
+  def __init__(self, Lgal, x, y=0.0, sculpt=-1):
-    "Support for Isothermal profile"
+    """Support for Isothermal profile
 
-    super(KSZ_Isothermal,self).__init__()
+    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__(sculpt)
     
     self.R_gal = 0.226 * x
     self.R_innergal *= y
@@ -71,7 +93,7 @@ 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)
@@ -84,7 +106,7 @@ class KSZ_Isothermal(KSZ_Profile):
     Q[cond] = ne.evaluate('rho0*2/(Mpc*r) * arctan(sqrt( (rGalaxy/r)**2 -1 ))',
                 local_dict=D)
 
-    return Q,np.where(r<rInner)[0]
+    return Q,[] #np.where(r<rInner)[0]
 
 
 # -----------------------------------------------------------------------------