Include all the stupid tools required for BORG analysis and plotting

python_sample/build_2lpt_ksz.py

@@ -0,0 +1,111 @@
+import healpy as hp
+import numpy as np
+import cosmotool as ct
+import argparse
+import h5py as h5
+from matplotlib import pyplot as plt
+
+Mpc=3.08e22
+rhoc = 1.8864883524081933e-26  # m^(-3)
+sigmaT = 6.6524e-29
+mp = 1.6726e-27
+lightspeed = 299792458.
+v_unit = 1e3  # Unit of 1 km/s
+T_cmb=2.725
+h = 0.71
+Y = 0.245  #The Helium abundance
+Omega_matter = 0.26
+Omega_baryon=0.0445
+
+G=6.67e-11
+MassSun=2e30
+frac_electron = 1.0 # Hmmmm
+frac_gas_galaxy = 0.14
+mu = 1/(1-0.5*Y)
+
+baryon_fraction = Omega_baryon / Omega_matter
+
+ksz_normalization = T_cmb*sigmaT*v_unit/(lightspeed*mu*mp) * baryon_fraction
+rho_mean_matter = Omega_matter * (3*(100e3/Mpc)**2/(8*np.pi*G)) 
+
+parser=argparse.ArgumentParser(description="Generate Skymaps from CIC maps")
+parser.add_argument('--boxsize', type=float, required=True)
+parser.add_argument('--Nside', type=int, default=128)
+parser.add_argument('--base_h5', type=str, required=True)
+parser.add_argument('--base_fig', type=str, required=True)
+parser.add_argument('--start', type=int, required=True)
+parser.add_argument('--end', type=int, required=True)
+parser.add_argument('--step', type=int, required=True)
+parser.add_argument('--minval', type=float, default=-0.5)
+parser.add_argument('--maxval', type=float, default=0.5)
+parser.add_argument('--depth_min', type=float, default=10)
+parser.add_argument('--depth_max', type=float, default=60)
+parser.add_argument('--iid', type=int, default=0)
+parser.add_argument('--ksz_map', type=str, required=True)
+args = parser.parse_args()
+
+L = args.boxsize
+Nside = args.Nside
+
+def build_sky_proj(density, dmax=60.,dmin=0,iid=0):
+
+  N = density.shape[0]
+  ix = (np.arange(N)-0.5)*L/N - 0.5 * L
+
+
+  dist2 = (ix[:,None,None]**2 + ix[None,:,None]**2 + ix[None,None,:]**2)
+
+  flux = density.transpose().astype(ct.DTYPE) # / dist2
+  dmax=N*dmax/L
+  dmin=N*dmin/L
+  if iid == 0:
+    shifter = np.array([0.5,0.5,0.5])
+  else:
+    shifter = np.array([0.,0.,0.])
+    
+  projsky1 = ct.spherical_projection(Nside, flux, dmin, dmax, integrator_id=iid, shifter=shifter, progress=1)
+
+  return projsky1*L/N
+
+def build_unit_vectors(N):
+  ii = np.arange(N)*L/N - 0.5*L
+  d = np.sqrt(ii[:,None,None]**2 + ii[None,:,None]**2 + ii[None,None,:]**2)
+  d[N/2,N/2,N/2] = 1
+  ux = ii[:,None,None] / d
+  uy = ii[None,:,None] / d
+  uz = ii[None,None,:] / d
+
+  return ux,uy,uz
+
+for i in xrange(args.start,args.end,args.step):
+  ff=plt.figure(1)
+  plt.clf()
+  with h5.File(args.base_h5 % i, mode="r") as f:
+    p = f['velocity'][:]
+    davg = np.average(np.average(np.average(f['density'][:],axis=0),axis=0),axis=0)
+    p /= davg # Now we have momentum scaled to the mean density
+
+  # Rescale by Omega_b / Omega_m
+  p = p.astype(np.float64)
+  print p.max(), p.min(), ksz_normalization, rho_mean_matter
+  p *= ksz_normalization*rho_mean_matter
+  p *= 1e6 # Put it in uK
+  p *= -1 # ksz has a minus
+  
+  ux,uy,uz = build_unit_vectors(p.shape[0])
+  pr = p[:,:,:,0] * ux + p[:,:,:,1] * uy + p[:,:,:,2] * uz 
+  print p.max(), p.min()
+  print pr.max()*Mpc, pr.min()*Mpc
+  
+  @ct.timeit_quiet
+  def run_proj():
+    return build_sky_proj(pr*Mpc, dmin=args.depth_min,dmax=args.depth_max,iid=args.iid)
+
+  proj = run_proj()
+
+  hp.write_map(args.ksz_map % i, proj)
+
+  hp.mollview(proj, fig=1, coord='CG', cmap=plt.cm.coolwarm, title='Sample %d' % i, min=args.minval, 
+              max=args.maxval)
+
+  ff.savefig(args.base_fig % i)

python_sample/build_2lpt_skymap.py

@@ -0,0 +1,68 @@
+import healpy as hp
+import numpy as np
+import cosmotool as ct
+import argparse
+import h5py as h5
+from matplotlib import pyplot as plt
+
+parser=argparse.ArgumentParser(description="Generate Skymaps from CIC maps")
+parser.add_argument('--boxsize', type=float, required=True)
+parser.add_argument('--Nside', type=int, default=128)
+parser.add_argument('--base_cic', type=str, required=True)
+parser.add_argument('--base_fig', type=str, required=True)
+parser.add_argument('--start', type=int, required=True)
+parser.add_argument('--end', type=int, required=True)
+parser.add_argument('--step', type=int, required=True)
+parser.add_argument('--minval', type=float, default=0)
+parser.add_argument('--maxval', type=float, default=60)
+parser.add_argument('--depth_min', type=float, default=10)
+parser.add_argument('--depth_max', type=float, default=60)
+parser.add_argument('--iid', type=int, default=0)
+args = parser.parse_args()
+
+INDATA="/nethome/lavaux/Copy/PlusSimulation/BORG/Input_Data/2m++.npy"
+tmpp = np.load(INDATA)
+
+L = args.boxsize
+Nside = args.Nside
+
+def build_sky_proj(density, dmax=60.,dmin=0,iid=0):
+
+  N = density.shape[0]
+  ix = (np.arange(N)-0.5)*L/N - 0.5 * L
+
+
+  dist2 = (ix[:,None,None]**2 + ix[None,:,None]**2 + ix[None,None,:]**2)
+
+  flux = density.transpose().astype(ct.DTYPE) # / dist2
+  dmax=N*dmax/L
+  dmin=N*dmin/L
+  if iid == 0:
+    shifter = np.array([0.5,0.5,0.5])
+  else:
+    shifter = np.array([0.,0.,0.])
+    
+  projsky1 = ct.spherical_projection(Nside, flux, dmin, dmax, integrator_id=iid, shifter=shifter)
+
+  return projsky1*L/N
+
+l,b = tmpp['gal_long'],tmpp['gal_lat']
+
+l = np.radians(l)
+b = np.pi/2 - np.radians(b)
+
+dcmb = tmpp['velcmb']/100.
+
+idx = np.where((dcmb>args.depth_min)*(dcmb<args.depth_max))
+
+
+for i in xrange(args.start,args.end,args.step):
+  ff=plt.figure(1)
+  plt.clf()
+  d = np.load(args.base_cic % i)
+  proj = build_sky_proj(1+d, dmin=args.depth_min,dmax=args.depth_max,iid=args.iid)
+
+  hp.mollview(proj, fig=1, coord='CG', cmap=plt.cm.coolwarm, title='Sample %d' % i, min=args.minval, max=args.maxval)
+  hp.projscatter(b[idx], l[idx], lw=0, color='g', s=2.0, alpha=0.7)
+
+  ff.savefig(args.base_fig % i)

python_sample/build_nbody_skymap.py

@@ -1 +0,0 @@
-/nethome/lavaux/wuala/WualaDrive/g_lavaux/PythonCode/BORG/build_nbody_skymap.py

python_sample/build_nbody_skymap.py

@@ -0,0 +1,58 @@
+import healpy as hp
+import numpy as np
+import cosmotool as ct
+import h5py as h5
+from matplotlib import pyplot as plt
+
+L=600.
+Nside=128
+
+INDATA="/nethome/lavaux/Copy/PlusSimulation/BORG/Input_Data/2m++.npy"
+tmpp = np.load(INDATA)
+
+def build_sky_proj(density, dmax=60.,dmin=0):
+
+  N = density.shape[0]
+  ix = (np.arange(N)-0.5)*L/N - 0.5 * L
+
+
+  dist2 = (ix[:,None,None]**2 + ix[None,:,None]**2 + ix[None,None,:]**2)
+
+  flux = density.transpose().astype(ct.DTYPE) # / dist2
+  dmax=N*dmax/L
+  dmin=N*dmin/L
+  projsky1 = ct.spherical_projection(Nside, flux, dmin, dmax, integrator_id=1)
+#  projsky0 = ct.spherical_projection(Nside, flux, 0, 52, integrator_id=0)
+
+  return projsky1*L/N#,projsky0
+
+l,b = tmpp['gal_long'],tmpp['gal_lat']
+
+l = np.radians(l)
+b = np.pi/2 - np.radians(b)
+
+dcmb = tmpp['velcmb']/100.
+
+idx = np.where((dcmb>10)*(dcmb<60))
+
+plt.figure(1)
+plt.clf()
+if True:
+  with h5.File("fields.h5", mode="r") as f:
+    d = f["density"][:].transpose()
+    d /= np.average(np.average(np.average(d,axis=0),axis=0),axis=0)
+    proj = build_sky_proj(d, dmin=10,dmax=60.)
+  proj0 = proj1 = proj
+else:
+  d = np.load("icgen/dcic0.npy")
+  proj0 = build_sky_proj(1+d, dmin=10,dmax=60.)
+  d = np.load("icgen/dcic1.npy")
+  proj1 = build_sky_proj(1+d, dmin=10,dmax=60.)
+
+hp.mollview(proj0, fig=1, coord='CG', max=60, cmap=plt.cm.coolwarm)
+hp.projscatter(b[idx], l[idx], lw=0, color='g', s=5.0, alpha=0.8)
+
+plt.figure(2)
+plt.clf()
+hp.mollview(proj1, fig=2, coord='CG', max=60, cmap=plt.cm.coolwarm)
+hp.projscatter(b[idx], l[idx], lw=0, color='g', s=5.0, alpha=0.8)

python_sample/gen_2lpt_asmooth.py

@@ -0,0 +1,47 @@
+import os
+import h5py as h5
+import numpy as np
+import cosmotool as ct
+import icgen as bic
+import icgen.cosmogrowth as cg
+import sys
+import argparse
+
+ADAPT_SMOOTH="/home/bergeron1NS/lavaux/Software/cosmotool/build/sample/simple3DFilter"
+cosmo={'omega_M_0':0.3175, 'h':0.6711}
+cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
+cosmo['omega_k_0'] = 0
+cosmo['omega_B_0']=0.049
+cosmo['SIGMA8']=0.8344
+cosmo['ns']=0.9624
+N0=256
+
+doSimulation=True
+
+astart=1.
+
+parser=argparse.ArgumentParser(description="Generate CIC density from 2LPT")
+parser.add_argument('--start', type=int, required=True)
+parser.add_argument('--end', type=int, required=True)
+parser.add_argument('--step', type=int, required=True)
+parser.add_argument('--base', type=str, required=True)
+parser.add_argument('--N', type=int, default=256)
+parser.add_argument('--output', type=str, default="fields_%d.h5")
+parser.add_argument('--supersample', type=int, default=1)
+args = parser.parse_args()
+
+
+
+for i in xrange(args.start, args.end, args.step):
+  print i
+
+  pos,vel,density,N,L,_ = bic.run_generation("%s/initial_density_%d.dat" % (args.base,i), 0.001, astart, 
+                                           cosmo, supersample=args.supersample, do_lpt2=True)
+
+  q = pos + vel
+
+  with h5.File("particles.h5", mode="w") as f:
+    f.create_dataset("particles", data=np.array(q).transpose())
+  
+  os.system(ADAPT_SMOOTH + " %s %lg %lg %d %lg %lg %lg" % ("particles.h5", 3000000, L, args.N, 0, 0, 0))
+  os.rename("fields.h5", args.output % i)

python_sample/gen_2lpt_density.py

@@ -0,0 +1,42 @@
+import numpy as np
+import cosmotool as ct
+import icgen as bic
+import icgen.cosmogrowth as cg
+import sys
+import argparse
+
+cosmo={'omega_M_0':0.3175, 'h':0.6711}
+cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
+cosmo['omega_k_0'] = 0
+cosmo['omega_B_0']=0.049
+cosmo['SIGMA8']=0.8344
+cosmo['ns']=0.9624
+N0=256
+
+doSimulation=True
+
+astart=1.
+
+parser=argparse.ArgumentParser(description="Generate CIC density from 2LPT")
+parser.add_argument('--start', type=int, required=True)
+parser.add_argument('--end', type=int, required=True)
+parser.add_argument('--step', type=int, required=True)
+parser.add_argument('--base', type=str, required=True)
+parser.add_argument('--N', type=int, default=256)
+parser.add_argument('--output', type=str, default="dcic_%d.npy")
+parser.add_argument('--supersample', type=int, default=1)
+args = parser.parse_args()
+
+
+
+for i in xrange(args.start, args.end, args.step):
+  print i
+#  pos,_,density,N,L,_ = bic.run_generation("/nethome/lavaux/remote/borg_2m++_128/initial_density_%d.dat" % i, 0.001, astart, cosmo, supersample=2, do_lpt2=True)
+  pos,_,density,N,L,_ = bic.run_generation("%s/initial_density_%d.dat" % (args.base,i), 0.001, astart, 
+                                           cosmo, supersample=args.supersample, do_lpt2=True)
+
+  dcic = ct.cicParticles(pos, L, args.N)
+  dcic /= np.average(np.average(np.average(dcic, axis=0), axis=0), axis=0)
+  dcic -= 1
+
+  np.save(args.output % i, dcic)

python_sample/icgen/test_whitify.py

@@ -0,0 +1,41 @@
+import numpy as np
+import cosmotool as ct
+import borgicgen as bic
+from matplotlib import pyplot as plt
+
+cosmo={'omega_M_0':0.3175, 'h':0.6711}
+cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
+cosmo['omega_k_0'] = 0
+cosmo['omega_B_0']=0.049
+cosmo['SIGMA8']=0.8344
+cosmo['ns']=0.9624
+
+zstart=50
+astart=1/(1.+zstart)
+halfPixelShift=False
+
+posx,vel,density,N,L,a_ic,cosmo = bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=1, shiftPixel=halfPixelShift, do_lpt2=False)
+
+w1 = bic.whitify(density, L, cosmo, supergenerate=1)
+w2 = bic.whitify(density, L, cosmo, supergenerate=2)
+
+N = w1.shape[0]
+Ns = w2.shape[0]
+
+w1_hat = np.fft.rfftn(w1)*(L/N)**3
+w2_hat = np.fft.rfftn(w2)*(L/Ns)**3
+
+P1, b1, dev1 = bic.bin_power(np.abs(w1_hat)**2, L, range=(0,3),bins=150,dev=True)
+P2, b2, dev2 = bic.bin_power(np.abs(w2_hat)**2, L, range=(0,3),bins=150,dev=True)
+
+fig = plt.figure(1)
+fig.clf()
+plt.fill_between(b1, P1*(1-dev1), P1*(1+dev1), label='Supergen=1', color='b')
+plt.fill_between(b2, P2*(1-dev2), P2*(1+dev2), label='Supergen=2', color='g', alpha=0.5)
+ax = plt.gca()
+ax.set_xscale('log')
+plt.ylim(0.5,1.5)
+plt.xlim(1e-2,4)
+plt.axhline(1.0, color='red', lw=4.0)
+plt.legend()
+plt.show()

python_sample/ramsesToArray.py

@@ -0,0 +1,14 @@
+import numpy as np
+import sys
+import h5py as h5
+import cosmotool as ct
+
+s = ct.loadRamsesAll(sys.argv[1], int(sys.argv[2]), doublePrecision=True, loadVelocity=True)
+
+
+q = [p for p in s.getPositions()]
+q += [p for p in s.getVelocities()]
+q = np.array(q)
+
+with h5.File("particles.h5", mode="w") as f:
+  f.create_dataset("particles", data=q.transpose())

python_sample/test_spheric_proj.py

@@ -1 +0,0 @@
-/nethome/lavaux/wuala/WualaDrive/g_lavaux/PythonCode/BORG/test_spheric_proj.py

python_sample/test_spheric_proj.py

@@ -0,0 +1,20 @@
+import cosmotool as ct
+import numpy as np
+import healpy as hp
+
+d = np.zeros((64,64,64), ct.DTYPE)
+
+d[32,32,32] = 1
+
+ii=np.arange(256)*64/256.-32
+xx = ii[:,None,None].repeat(256,axis=1).repeat(256,axis=2).reshape(256**3)
+yy = ii[None,:,None].repeat(256,axis=0).repeat(256,axis=2).reshape(256**3)
+zz = ii[None,None,:].repeat(256,axis=0).repeat(256,axis=1).reshape(256**3)
+
+d_high = ct.interp3d(xx, yy, zz, d, 64, periodic=True)
+d_high = d_high.reshape((256,256,256))
+
+proj0 = ct.spherical_projection(64, d, 0, 20, integrator_id=0, shifter=np.array([0.5,0.5,0.5]))
+proj1 = ct.spherical_projection(64, d, 0, 20, integrator_id=1)
+
+proj0_high = ct.spherical_projection(256, d_high, 0, 30, integrator_id=0, shifter=np.array([0.5,0.5,0.5]))