Merge branch 'dev' into joss

README.md

@@ -1,6 +1,6 @@
 # `DeepSZSim`
 
-Code for producing fast simulations of the SZ effect for galaxy halos of varying z, $M_{200}$, based on average thermal pressure profile fits from [Battaglia et al. 2012](https://ui.adsabs.harvard.edu/abs/2012ApJ...758...75B/abstract). Simulated submaps can include tSZ signal from these halos, simulated CMB, instrument beam convolution and white noise.
+Code for producing fast simulations of the SZ effect for galaxy halos of varying redshift and mass, based on average thermal pressure profile fits from [Battaglia et al. 2012](https://ui.adsabs.harvard.edu/abs/2012ApJ...758...75B/abstract). Simulated submaps can include tSZ signal from these halos, simulated CMB, instrument beam convolution and white noise.
 
 ## Code Overview
 
@@ -11,19 +11,19 @@ The CMB simulations are handled by [DeepCMBSim](https://github.com/deepskies
 
 ### Installation 
 
-We provide an environment specification file for `conda` or `mamba` users at `environment.yml`. With `conda`, an environment is created by `conda env create -f environment.yml`. With `micromamba` the `env` is omitted and a new environment is instead created with `micromamba create -f environment.yml`.
+We provide an environment specification file for `conda` or `mamba` users at `environment.yml`, which will produce a new virtual environment called `szsims` with appropriate versions of major python packages. With `conda`, this environment is created by `conda env create -f environment.yml`. With `micromamba` the `env` is omitted and a new environment is instead created with `micromamba create -f environment.yml`.
 
-The simulated CMB signal relies on `camb` and utilities for saving rely on `h5py`.
+The simulated CMB signal relies on `camb` and `pixell`, cosmology relies on `colossus`, and utilities for saving rely on `h5py`. These are specified in the `pyproject.toml` file.
 
-From the top-level directory, you can do `pip install .`
+From the top-level directory, you can do `pip install .` to install the package.
 
 ### Usage
 
 The usage of this code is documented in `notebooks/demo_simulation.ipynb`. A detailed walkthrough of the functions available in this code is in `notebooks/demo_full_pipeline.ipynb`.
 
 A full list of potential inputs is documented in `settings/config.yaml` and you can edit `settings/inputdata.yaml` to reflect your desired simulation settings.  
 
-`dm_halo_dist.py` generates a z, $M_{200}$ array. The functions in `make_sz_cluster.py` create pressure profiles, Compton-y, and SZ signal maps from these halos of various z, $M_{200}$ and produce the final simulated submaps. These submaps contain simulated CMB and simple instrument beam convolution from `simtools.py` and white noise from `noise.py`. Plotting tools are provided in `visualization.py`.
+`dm_halo_dist.py` generates an array of mass and redshift. The functions in `make_sz_cluster.py` create pressure profiles, Compton-y, and SZ signal maps for a halo of a given mass and redshift, and produces the final simulated submaps. These submaps contain simulated CMB and simple instrument beam convolution from `simtools.py` and white noise from `noise.py`. Plotting tools are provided in `visualization.py`. Simulations of a large suite of clusters can be achieved easily with `simclusters.py`.
 
 ### Example
 
@@ -43,7 +43,7 @@ To access the clusters in this set, you can refer to the cluster ID, which itsel
 
 ## Citation
 
-If you use this code in your research, please cite this GitHub repo. Please also make use of the citation instructions for `camb` provided [here](https://camb.info).
+If you use this code in your research, please cite this GitHub repo and our JOSS paper. Please also make use of the citation instructions for `camb` provided [here](https://camb.info).
 
 ## Contributing
 

deepszsim/dm_halo_dist.py

@@ -3,8 +3,6 @@
 """
 
 import numpy as np
-from colossus.cosmology import cosmology
-from colossus.halo import mass_adv
 
 def flatdist_halo(zmin, zmax, m200min_SM, m200max_SM, size, seed=None):
     '''

deepszsim/filters.py

@@ -4,16 +4,18 @@
 
 import numpy as np
 
-def get_tSZ_signal_aperture_photometry(dT_map, radmax_arcmin, 
-                                           fmax_arcmin=np.sqrt(2)):
+def get_tSZ_signal_aperture_photometry(dT_map, radmax_arcmin, pixel_scale, 
+                                       fmax=np.sqrt(2)):
     """
     Parameters:
     ----------
     dT_map: array to represent the map in uK
     radmax_arcmin: float
         the radius of the inner aperture, in arcmin
-    fmax_arcmin: float
-        fmax+radmax is the radius of the outer aperture, in arcmin
+    pixel_scale: float
+        How many arcmin per pixel for the current settings
+    fmax: float
+        fmax * radmax_arcmin is the radius of the outer aperture, in arcmin
 
     Returns:
     -------
@@ -25,16 +27,14 @@ def get_tSZ_signal_aperture_photometry(dT_map, radmax_arcmin,
         thermal SZ effect signal
     """
 
-    center = np.array(dT_map.shape) / 2
+    radmax_pixels = radmax_arcmin / pixel_scale
+    radius_out_pixels = radmax_pixels * fmax
+
+    center = np.array(dT_map.shape) // 2
     x, y = np.indices(dT_map.shape)
     r = np.sqrt((x - center[0])**2 + (y - center[1])**2)
-
-    #define outer radius
-    radius_out=radmax_arcmin * fmax_arcmin 
-    #average in inner radius
-    disk_mean = dT_map[r < radmax_arcmin].mean() 
-    #average in outer radius
-    ring_mean = dT_map[(r >= radmax_arcmin) & (r < radius_out)].mean()
+
+    disk_mean = dT_map[r < radmax_pixels].mean()
+    ring_mean = dT_map[(r >= radmax_pixels) & (r < radius_out_pixels)].mean()
     tSZ = disk_mean - ring_mean
-
     return disk_mean, ring_mean, tSZ

deepszsim/make_sz_cluster.py

@@ -228,8 +228,7 @@ def Pth_Battaglia2012(radius_mpc, M200_SM, redshift_z, load_vars_dict = None,
     return _Pth_Battaglia2012(P0, radius_mpc, R200_Mpc, alpha, beta, gamma, xc)
 
 
-def Pe_to_y(profile, radii_mpc, M200_SM, redshift_z, load_vars_dict, alpha = 1.0, gamma = -0.3, R200_Mpc = None,
-            Rmaxy = None):
+def Pe_to_y(profile, radii_mpc, M200_SM, redshift_z, load_vars_dict, alpha = 1.0, gamma = -0.3, R200_Mpc = None):
     '''
     Converts from an electron pressure profile to a compton-y profile,
     integrates over line of sight from -1 to 1 Mpc relative to center.
@@ -264,33 +263,19 @@ def Pe_to_y(profile, radii_mpc, M200_SM, redshift_z, load_vars_dict, alpha = 1.0
     if R200_Mpc is None:
         R200_Mpc = get_r200_angsize_and_c200(M200_SM, redshift_z, load_vars_dict)[1]
     radii_mpc = (radii_mpc * u.Mpc).value
-    if Rmaxy is None:
-        rmax = radii_mpc.max()
-    elif '200' in Rmaxy:
-        rmax = R200_Mpc
-    else:
-        print('please specify a valid `Rmaxy`')
-        return None
     if profile != "Battaglia2012":
         print("only implementing `Battaglia2012` for profile")
     profile = Pth_Battaglia2012
     pressure_integ = np.empty_like(radii_mpc)
     P200_kevcm3 = P200_Battaglia2012(M200_SM, redshift_z, load_vars_dict, R200_Mpc = R200_Mpc).value
 
-    # integral = np.trapz(np.array([profile(np.sqrt(np.linspace(0, np.sqrt(radii_mpc.max()**2. - rv**2.)+1.,
-    #                                                               1000)**2 +
-    #                                      rv**2), M200_SM, redshift_z, load_vars_dict = None, alpha = alpha,
-    #                       gamma = gamma, R200_Mpc = r200) for rv in radii_mpc]), np.array([np.linspace(0,
-    #                                                                                             np.sqrt(radii_mpc.max(
-    # )**2. - rv**2.)+1., 1000) for rv in radii_mpc]))
-    # y_pro = integral * P200_kevcm3 * keVcm3_to_Jm3 * Thomson_scale * \
-    #         thermal_to_electron_pressure * 2*Mpc_to_m
-    # return y_pro
+    rmax = radii_mpc.max()
+
     for i, radius in enumerate(radii_mpc):
         # Multiply profile by P200 specifically for Battaglia 2012 profile,
         # since it returns Pth/P200 instead of Pth
         rv = radius
-        if (rmax == R200_Mpc) and (rv >= R200_Mpc):
+        if (rv >= R200_Mpc):
             pressure_integ[i] = 0
         else:
             l_mpc = np.linspace(0, np.sqrt(rmax**2. - rv**2.) + 1., 1000)  # Get line of sight axis
@@ -303,7 +288,7 @@ def Pe_to_y(profile, radii_mpc, M200_SM, redshift_z, load_vars_dict, alpha = 1.0
 
 
 def _make_y_submap(profile, M200_SM, redshift_z, load_vars_dict, image_size_pixels, pixel_size_arcmin, alpha = 1.0,
-                   gamma = -0.3, R200_Mpc = None, Rmaxy = None):
+                   gamma = -0.3, R200_Mpc = None):
     '''
     Converts from an electron pressure profile to a compton-y profile,
     integrates over line of sight from -1 to 1 Mpc relative to center.
@@ -324,6 +309,10 @@ def _make_y_submap(profile, M200_SM, redshift_z, load_vars_dict, image_size_pixe
         size of final submap in number of pixels
     pixel_size_arcmin: float
         size of each pixel in arcmin
+    alpha: float
+        variable fixed by Battaglia et al 2012 to 1.0
+    gamma: float
+        variable fixed by Battaglia et al 2012 to -0.3
     R200_Mpc: None or float
         if None, will calculate the radius that corresponds to the mass M200, the redshift redshift_z,
         and the cosmology contained in load_vars_dict
@@ -341,8 +330,7 @@ def _make_y_submap(profile, M200_SM, redshift_z, load_vars_dict, image_size_pixe
     mindist = utils.arcmin_to_Mpc(pixel_size_arcmin*0.1, redshift_z, load_vars_dict['cosmo'])
     R = np.maximum(mindist, np.sqrt(X[:, None]**2 + X[None, :]**2).flatten())
 
-    cy = Pe_to_y(profile, R, M200_SM, redshift_z, load_vars_dict, alpha = alpha, gamma = gamma, R200_Mpc = R200_Mpc,
-                 Rmaxy = Rmaxy)  #
+    cy = Pe_to_y(profile, R, M200_SM, redshift_z, load_vars_dict, alpha = alpha, gamma = gamma, R200_Mpc = R200_Mpc)  #
     # evaluate compton-y for each
     # neccesary radius
 
@@ -366,7 +354,7 @@ def _make_y_submap(profile, M200_SM, redshift_z, load_vars_dict, image_size_pixe
 
 def generate_y_submap(M200_SM, redshift_z, profile = "Battaglia2012",
                       image_size_pixels = None, pixel_size_arcmin = None, load_vars_dict = None, alpha = 1.0, gamma = -0.3,
-                      R200_Mpc = None, Rmaxy = None):
+                      R200_Mpc = None):
     '''
     Converts from an electron pressure profile to a compton-y profile,
     integrates over line of sight from -1 to 1 Mpc relative to center.
@@ -387,6 +375,10 @@ def generate_y_submap(M200_SM, redshift_z, profile = "Battaglia2012",
     load_vars_dict: dict
         result of running the load_vars() function, which includes a dictionary of cosmological and experimental
         parameters
+    alpha: float
+        variable fixed by Battaglia et al 2012 to 1.0
+    gamma: float
+        variable fixed by Battaglia et al 2012 to -0.3
     R200_Mpc: None or float
         if None, will calculate the radius that corresponds to the mass M200, the redshift redshift_z,
         and the cosmology contained in load_vars_dict
@@ -406,11 +398,11 @@ def generate_y_submap(M200_SM, redshift_z, profile = "Battaglia2012",
 
     y_map = _make_y_submap(profile, M200_SM, redshift_z, load_vars_dict,
                            image_size_pixels, pixel_size_arcmin,
-                           alpha = alpha, gamma = gamma, R200_Mpc = R200_Mpc, Rmaxy = Rmaxy)
+                           alpha = alpha, gamma = gamma, R200_Mpc = R200_Mpc)
 
     return y_map
 
-def get_r200_angsize_and_c200(M200_SM, redshift_z, load_vars_dict):
+def get_r200_angsize_and_c200(M200_SM, redshift_z, load_vars_dict, angsize_density = None):
     '''
     Parameters:
     ----------
@@ -421,6 +413,9 @@ def get_r200_angsize_and_c200(M200_SM, redshift_z, load_vars_dict):
     load_vars_dict: dict
         must contain 'cosmo' (a FlatLambaCDM instance describing the background cosmology),
         'sigma8' (float, around 0.8), and 'ns' (float, around 0.96)
+    angsize_density: None or str
+        density measure at which to calculate the angular size, if desired. If `None`, will not
+        calculate an angular size. Otherwise, use a valid choice as specified in `colossus.halo.mass_adv`
 
     Returns:
     -------
@@ -439,9 +434,19 @@ def get_r200_angsize_and_c200(M200_SM, redshift_z, load_vars_dict):
     cosmology.addCosmology('myCosmo', **params)
     cosmo_colossus = cosmology.setCosmology('myCosmo')
 
-    M200_SM, R200_Mpc, c200 = mass_adv.changeMassDefinitionCModel(M200_SM * cosmo.h,
+    M200_SM, R200_kpc, c200 = mass_adv.changeMassDefinitionCModel(M200_SM * cosmo.h,
                                                                   redshift_z, '200c', '200c', c_model = 'ishiyama21')
     M200_SM /= cosmo.h  # From M_solar/h to M_solar
-    R200_Mpc = R200_Mpc / cosmo.h / 1000  # From kpc/h to Mpc
-    angsize_arcmin = R200_Mpc*1000/60/cosmo_colossus.kpcPerArcsec(redshift_z)
+    R200_Mpc = R200_kpc / cosmo.h / 1000  # From kpc/h to Mpc
+
+    if angsize_density is not None:
+        if angsize_density != '200c':
+            _, Rd_kpc, _ = mass_adv.changeMassDefinitionCModel(M200_SM * cosmo.h,
+                                                               redshift_z, '200c', angsize_density,
+                                                               c_model = 'ishiyama21')
+            angsize_arcmin = Rd_kpc / cosmo.h / 60 / cosmo_colossus.kpcPerArcsec(redshift_z)
+        else:
+            angsize_arcmin = R200_Mpc * 1000 / 60 / cosmo_colossus.kpcPerArcsec(redshift_z)
+    else:
+        angsize_arcmin = None
     return M200_SM, R200_Mpc, angsize_arcmin, c200 # now returns M200, R200, angsize in arcmin, c200

deepszsim/noise.py

@@ -9,8 +9,12 @@ def generate_noise_map(image_size, noise_level, pix_size):
     Generates a white noise map based on the noise level and beam size.
 
     Args:
-        image_size (int): Size of the noise map (N x N).
-        noise_level (float): Noise level of the survey.
+        image_size: int
+            Size of the noise map (N x N).
+        noise_level: float
+            Noise level of the survey.
+        pix_size: int
+            size of pixels in arcminutes
 
     Returns:
         ndarray: Noise map.

deepszsim/simclusters.py

@@ -10,7 +10,7 @@
 
 class simulate_clusters:
     def __init__(self, M200 = None, redshift_z = None, num_halos = None, halo_params_dict = None,
-                 R200_Mpc = None, Rmaxy = None, profile = "Battaglia2012",
+                 R200_Mpc = None, profile = "Battaglia2012",
                  image_size_pixels = None, image_size_arcmin = None, pixel_size_arcmin = None,
                  alpha = 1.0, gamma = -0.3,
                  load_vars_yaml = os.path.join(os.path.dirname(__file__), 'Settings', 'inputdata.yaml'),
@@ -19,9 +19,9 @@ def __init__(self, M200 = None, redshift_z = None, num_halos = None, halo_params
         """
         Parameters
         ----------
-        M200_dist: float or array-like of float
+        M200: float or array-like of float
             the mass contained within R200 in solar masses (same length as z_dist)
-        z_dist: float or array-like of float
+        redshift_z: float or array-like of float
             the redshift of the cluster (unitless) (same length as M200_dist)
         num_halos: None or int
             number of halos to simulate if none supplied
@@ -47,6 +47,8 @@ def __init__(self, M200 = None, redshift_z = None, num_halos = None, halo_params
             path to yaml file with params
         seed: None or int
             random seed value to sample with
+        tqverb: bool
+            whether or not to display tqdm progress bar while making T maps
         """
 
         if (M200 is not None) and (redshift_z is not None):
@@ -93,25 +95,23 @@ def __init__(self, M200 = None, redshift_z = None, num_halos = None, halo_params
         if R200_Mpc is not None:
             self.R200_Mpc = R200_Mpc
         else:
-            self.R200_Mpc, self.angsize_arcmin = np.array(
-                [make_sz_cluster.get_r200_angsize_and_c200(self.M200[i], self.redshift_z[i], self.vars)[1:-1]
+            self.R200_Mpc, self.angsize500_arcmin = np.array(
+                [make_sz_cluster.get_r200_angsize_and_c200(self.M200[i], self.redshift_z[i], self.vars,
+                                                           angsize_density = '500c')[1:3]
                  for i in range(self._size)]).T
 
-        self.Rmaxy = Rmaxy
-
         self.id_list = [
             str(self.M200[i])[:5] + str(self.redshift_z[i] * 100)[:2] + str(self._rng.integers(10**6)).zfill(6)
             for i in range(self._size)]
         self.clusters.update(zip(self.id_list, [{"params": {'M200': self.M200[i], 'redshift_z': self.redshift_z[i],
                                                             'R200': self.R200_Mpc[i],
-                                                            'angsize_arcmin': self.angsize_arcmin[i],
+                                                            'angsize500_arcmin': self.angsize500_arcmin[i],
                                                             'image_size_pixels': self.image_size_pixels}} for
                                                 i in range(
                 self._size)]))
 
     def get_y_maps(self):
         """
-
         Returns
         -------
         np.ndarray(float)
@@ -124,14 +124,12 @@ def get_y_maps(self):
             self.y_maps = np.array([make_sz_cluster.generate_y_submap(self.M200[i],
                                                       self.redshift_z[i],
                                                       R200_Mpc = self.R200_Mpc[i],
-                                                      Rmaxy = self.Rmaxy,
                                                       load_vars_dict = self.vars)
                                     for i in tqdm(range(self._size), disable = (not self.tqverb))])
             return self.y_maps
 
     def get_dT_maps(self):
         """
-
         Returns
         -------
         np.ndarray(float)
@@ -152,6 +150,8 @@ def get_T_maps(self, add_CMB = True, returnval = False):
         ----------
         add_CMB: bool
             whether or not to include the CMB contribution to the final map
+        returnval: bool
+            whether or not to return the T maps themselves or simply update internal attribute
 
         Returns
         -------

environment.yml

@@ -1,5 +1,6 @@
 name: szsims
 dependencies:
+  - python < 3.12
   - numpy
   - pandas
   - astropy

notebooks/demo_simulation.ipynb

@@ -266,9 +266,9 @@
    "provenance": []
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Any Display Name",
    "language": "python",
-   "name": "python3"
+   "name": "your_env_name"
   },
   "language_info": {
    "codemirror_mode": {
@@ -280,7 +280,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.11.0"
   }
  },
  "nbformat": 4,