Overview of the code¶
[1]:
%pylab inline
from pmesh.pm import ParticleMesh
from gaussianCR.construct import *
from gaussianCR.cosmo import *
Populating the interactive namespace from numpy and matplotlib
[2]:
np.set_printoptions(precision=2,linewidth=150,suppress=True)
initialize cosmology
First, we need to initialize a Cosmos object that collects some variables and functions dependent on the cosmological parameters.
One option (recommendated) to initialize a Cosmos object is using nbodykit.cosmology.
[3]:
import nbodykit.cosmology as nbcosmos
wmap9 = Cosmos(FLRW=True,obj=nbcosmos.WMAP9)
If set FLRW==True, we use nbodykit.cosmology.LinearPower to compute the linear power spectrum for the chosen cosmology, using the analytic Eisenstein & Hu approximation. Note that obj can also be an instance from astropy.cosmology.FLRW.
Alternatively, one can also intialize Cosmos by setting flag FLRW = Flase, and manually input the cosmological parameters and linear matter power spectrum (at z=0). Linear power spectrum can be the output generated from CAMB, CLASS or some alternative cosmology calculater. Note that it should be consistent with the input cosmological parameters.
e.g. Below we initialize Cosmos with fuzzy dark matter power spectrum generated by axionCAMB.
[4]:
data = np.loadtxt('/home/yueying/source/gaussianCR/examples/ps/fdm_m2p5_z0_matterpower.dat')
pk_func = interpolate.interp1d(data[:,0],data[:,1],fill_value='extrapolate')
mycosmos = Cosmos(FLRW=False,H0=69.3,Om0=0.286,Ob0=0.0463,Pk_lin=pk_func)
initialze gsCR
gsCR(cosmology, Lbox=20, Nmesh=128, RG=0.9, xpk=[0, 0, 0], CONS=[‘full’])
Here we initialize gsCR object with:
- wmap9 : Cosmos object we build above
- Lbox : Boxsize of the initial condition in Mpc/h
- Nmesh : the grid size to represent the linear density field, in shape of (Nmesh,Nmesh,Nmesh)
- RG : The size of the Gaussian kernel we use to impose the constraints, in unit of Mpc/h
- xpk : Position to impose the constraint, in unit of Mpc/h. Default is (0, 0, 0), we can change it later.
- CONS : The flags to turn on the specific constraints, details see below.
The available options of CONS in (‘full’,’f0’,’f1’,’f2’,’vx’,’vy’,’vz’,’TG’):
- full : enable all the 18 constraints at position xpk
- f0 : constrain the height of the density peak (zeroth order of fG field)
- f1 : constrain the three 1st order derivatives of fG field at xpk,
- f2 : constrain the the compactness,ellipticity and orientation of the peak (2nd order derivatives of fG field)
- vx,vy,vz : constrain the three peculiar velocities of fG field at xpk
- TG : constrain the tidal field of fG field at xpk
For example, below we initialize a gsCR object fg. We can set xpk and CONS later.
[5]:
fg = gsCR(wmap9,Lbox=20,Nmesh=128,RG=0.9)
print (fg)
This is a gsCR object:
Lbox = 20.0 Mpc/h
Nmesh = 128
RG = 0.9 Mpc/h
Sigma0_RG = 1.82, Sigma2_RG = 1.78
xpk = [0, 0, 0]
CONS = ['full']
xij_tensor_inv =
None
gsCR.build_Xij_inv_matrix()
Sometimes we only want to impose a subset of the full 18 peak constraints. For example, we only want to manipulate the height and shape of the density peak, but leave the tidal field and peculiar velocity free. Then we should set CONS to be [‘f0’,’f2’] (which corresponds to H0, H5 - H10) and build the corresponding \(\xi_{ij}^{-1}\) matrix. gsCR.build_Xij_inv_matrix() method would initialize the gsCR.xij_tensor_inv covariance matrix constructed by the kernels corresponding to
CONS flag. This function need to be called whenever we set or change the CONS flags.
Below we set the CONS = [‘f0’,’f2’]. The corresponding \(\xi_{ij}^{-1}\) matrix looks like:
[6]:
fg.CONS = ['f0','f2']
fg.build_Xij_inv_matrix()
print (fg.xij_tensor_inv)
[[ 0.63 0.38 0.38 0.38 -0. 0. 0. ]
[ 0.38 2.12 -0.25 -0.25 0. 0. -0. ]
[ 0.38 -0.25 2.12 -0.25 -0. 0. -0. ]
[ 0.38 -0.25 -0.25 2.12 0. -0. -0. ]
[-0. 0. -0. 0. 4.75 -0. 0. ]
[ 0. 0. 0. -0. -0. 4.75 -0. ]
[ 0. -0. -0. -0. 0. -0. 4.75]]
Here is the full \(\xi_{ij}^{-1}\) matrix with the 18 constraints:
[7]:
fg.CONS = ['full']
fg.build_Xij_inv_matrix()
print (fg.xij_tensor_inv)
[[ 0.63 0. 0. 0. 0.38 0.38 0.37 -0. 0. 0. 0. 0. 0. -0. -0. 0. -0. -0. ]
[ 0. 2.49 0. 0. 0. 0. 0. 0. 0. 0. -0.02 0. -0. 0. 0. 0. 0. 0. ]
[ 0. 0. 2.49 0. 0. 0. 0. 0. 0. 0. 0. -0.02 0. 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 2.49 0. 0. 0. 0. 0. 0. -0. 0. -0.02 0. 0. 0. 0. 0. ]
[ 0.38 0. 0. 0. 2.99 -0.68 -0.68 0. 0. -0. 0. 0. 0. -0.07 0. -0. -0. 0. ]
[ 0.38 0. 0. 0. -0.68 2.99 -0.68 -0. 0. 0. 0. 0. 0. 0. -0.07 0. -0. -0. ]
[ 0.37 0. 0. 0. -0.68 -0.68 2.99 0. -0. -0. 0. 0. 0. 0.06 0.06 -0. 0. 0. ]
[-0. 0. 0. 0. 0. -0. 0. 9.03 -0. -0. 0. 0. 0. -0. 0. -0.22 0. 0. ]
[ 0. 0. 0. 0. 0. 0. -0. -0. 9.03 0. 0. 0. 0. -0. -0. 0. -0.22 -0. ]
[ 0. 0. 0. 0. -0. 0. -0. -0. 0. 9.03 0. 0. 0. -0. -0. 0. -0. -0.22]
[-0. -0.02 -0. -0. -0. -0. -0. -0. -0. -0. 0. -0. 0. -0. -0. -0. -0. -0. ]
[-0. -0. -0.02 -0. -0. -0. -0. -0. -0. -0. -0. 0. -0. -0. -0. -0. -0. -0. ]
[-0. -0. -0. -0.02 -0. -0. -0. -0. -0. -0. 0. -0. 0. -0. -0. -0. -0. -0. ]
[-0. -0. -0. -0. -0.07 0. 0.06 -0. -0. -0. -0. -0. -0. 0.01 0. 0. 0. 0. ]
[-0. -0. -0. -0. 0. -0.07 0.06 0. -0. -0. -0. -0. -0. 0. 0.01 -0. 0. 0. ]
[ 0. -0. -0. -0. -0. 0. -0. -0.22 0. 0. -0. -0. -0. 0. -0. 0.01 -0. -0. ]
[-0. -0. -0. -0. -0. -0. 0. 0. -0.22 -0. -0. -0. -0. 0. 0. -0. 0.01 0. ]
[-0. -0. -0. -0. 0. -0. 0. 0. -0. -0.22 -0. -0. -0. 0. 0. -0. 0. 0.01]]
gsCR.set_c_values()
Apart from the \(\xi_{ij}^{-1}\) matrix, we also need to set {\(c_i\)} corresponding to the constraint kernels. gsCR.set_c_values() transfer the peak features to the corresponding {\(c_i\)} sets.
arguments:
- nu : peak height, in unit of \(\sigma_0(R_G)\)
- xd : peak compactness, in unit of \(\sigma_2(R_G)\)
- a12sq : axial ratio \((a_1/a_2)^2\)
- a13sq : axial ratio \((a_1/a_3)^2\)
- a1,b1,p1 : Euler angle to transform to principal axis of mass ellipsoid
- vx,vy,vz : Peculiar velocity of the peak in unit km/s
- epsilon : Shear magnitude in unit km/s/Mpc
- omega : Shear angle to distribute the shear magnitude between three axes, \([\pi,2\pi]\)
- a2,b2,p2 : Euler angle to transform to principal axis of tidal tensor
- silent : flag to print the relevant peak parameters
Note that we only need to set arguments corresponding to CONS flag, other will be ignored
gsCR.Ensemble_field()
With \(\xi_{ij}^{-1}\) matrix and \(c_i\) values ready, we can build the Ensemble mean field with \(\bar{f}(\mathbf x) = <f(\mathbf x)|\Gamma> = \xi_i(\mathbf x) \xi^{-1}_{ij} c_j\):
gsCR.Ensemble_field() takes the argument of {\(c_i\)} which is the c_value obtained from gsCR.set_c_values(), and return the density field of the ensemble mean field.
See next section for more details and examples.
gsCR.read_out_c18()
This is used to obtain the original {\(c_i\)} values of the unconstrained density field at position rpos.
argument: - dx_field : with shape (Ng,Ng,Ng), the density contrast field, the boxsize should match self.attrs[‘Lbox’] - rpos : the position to read out c values, if not set, then use self.xpk
return: - {\({c_i}\)} sets : with i=1,…18. - peak_data : structured array that converts the {\(c_i\)} sets to the peak parameters.