Merge pull request #237 from deepskies/dev

cleaning up small things
deepskies · Sep 28, 2024 · d0aa047 · d0aa047
2 parents d3a871c + 9888652
commit d0aa047
Show file tree

Hide file tree

Showing 28 changed files with 1,346 additions and 145 deletions.
diff --git a/.github/workflows/paper.yml → .github/workflows/draft-pdf.yml b/.github/workflows/paper.yml → .github/workflows/draft-pdf.yml
@@ -1,23 +1,28 @@
-on: [push]
+name: Draft PDF
+on:
+  push:
+    paths:
+      - paper/**
+      - .github/workflows/draft-pdf.yml
 
 jobs:
   paper:
     runs-on: ubuntu-latest
     name: Paper Draft
     steps:
       - name: Checkout
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
       - name: Build draft PDF
         uses: openjournals/openjournals-draft-action@master
         with:
           journal: joss
           # This should be the path to the paper within your repo.
           paper-path: paper/paper.md
       - name: Upload
-        uses: actions/upload-artifact@v1
+        uses: actions/upload-artifact@v3
         with:
           name: paper
           # This is the output path where Pandoc will write the compiled
           # PDF. Note, this should be the same directory as the input
           # paper.md
-          path: paper/paper.pdf
+          path: paper/paper.pdf
diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
@@ -9,18 +9,17 @@ jobs:
       max-parallel: 5
 
     steps:
-    - uses: actions/checkout@v3
     - name: Set up Python 3.10
       uses: actions/setup-python@v3
       with:
         python-version: '3.10'
-    - name: Add conda to system path
-      run: |
-        # $CONDA is an environment variable pointing to the root of the miniconda directory
-        echo $CONDA/bin >> $GITHUB_PATH
+    - uses: actions/checkout@v4
+    - uses: conda-incubator/setup-miniconda@v3
+      with:
+        activate-environment: base-test
     - name: Install dependencies
       run: |
-        conda env update --file=environment.yml --name=base
+        # conda env update --file=../environment.yml --name=base-test
         python -m pip install colossus pixell camb
         python -m pip install .
     - name: Test with pytest

diff --git a/.gitignore b/.gitignore
@@ -2,6 +2,8 @@
 __pycache__/
 *.py[cod]
 *$py.class
+**/*.tfevents.*
+.DS_Store
 
 # C extensions
 *.so

diff --git a/README.md b/README.md
@@ -1,26 +1,49 @@
 # `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.
 
-## Installation 
+## Code Overview
 
-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`.
+The code is structured as depicted here: ![DeepSZSim workflow](paper/figures/DeepSZSim_Workflow.png)
+The CMB simulations are handled by [DeepCMBSim](https://www.github.com/deepskies/deepcmbsim), based on [CAMB](https://camb.info), and further by [pixell](https://github.com/simonsobs/pixell). The SZ cluster simluations are done in `make_sz_cluster.py` and instrumental effects are added in `filters.py` and `noise.py`.  
 
-The simulated CMB signal relies on `camb` and utilities for saving rely on `h5py`.
+## Quickstart
 
-From the top-level directory, you can do `pip install .`
+### Installation 
 
-## Usage
+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 `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 .` 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
+
+Let's say you wanted to produce 100 mock halos distributed across the redshift range 0.2<z<0.4 and with masses in the range 1e14<M200<1e15. To generate these halos and produce their simulated maps with SZ signal (along with CMB signal and noise parameters as specified in `Settings/inputdata.yaml`) you would call
+```commandline
+import deepszsim as dsz
+tc0 = dsz.simulate_clusters(halo_params_dict={
+                            'zmin':0.2, 'zmax':0.5,
+                            'm200min_SM':1e14, 'm200max_SM':1e15
+                            },
+                            num_halos=100)
+tc0.get_T_maps()
+```
+The clusters and their maps are now in a dictionary which is in a `clusters` attribute of the class instance `tc0`.
+
+To access the clusters in this set, you can refer to the cluster ID, which itself is obtained from the first five digits of the cluster mass and two digits of the cluster redshift, followed by six random digits. For example, to access a dictionary of the maps and the parameters describing the eleventh cluster, you would do `tc0.clusters[tc0.id_list[11]]`. Alternately, to get the ''final'' temperature map (with noise) for the eleventh cluster, we also provide a convenience function: `tc0.ith_T_map(11)` is the same as `tc0.clusters[tc0.id_list[11]]['maps']['final_map']`.
 
 ## 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
 

diff --git a/deepszsim/Settings/config.yaml → deepszsim/Settings/config_everything.yaml b/deepszsim/Settings/config.yaml → deepszsim/Settings/config_everything.yaml
@@ -1,16 +1,3 @@
-#Reformat to accomodate decisions on astropy constants and cosmology param input 
-#UNIVERSAL_CONSTANTS:
-#  planck_const: 6.626e-34
-#  boltzman_const: 1.38e-23
-#  G: 6.6743e-11        #m^3/Kg/s^2
-#  m_sun: 1.98847e30       #Kg
-#  Thomson_sec: 6.65246e-29       #m^2
-#  m_electron: 9.11e-31        #Kg
-#  c: 299792458     #m/s
-#  Mpc_to_m: 3.09e22
-#  kevcm_to_jm: 1.6e-16 * 1e6
-#  j_to_kev: 6.242e15
-
 COSMOLOGY:
   Planck2018: #https://arxiv.org/abs/1807.06209:
     flat: True 

diff --git a/deepszsim/Settings/inputdata.yaml → deepszsim/Settings/config_simACTDR5.yaml b/deepszsim/Settings/inputdata.yaml → deepszsim/Settings/config_simACTDR5.yaml
diff --git a/deepszsim/dm_halo_dist.py b/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):
     '''
@@ -23,6 +21,8 @@ def flatdist_halo(zmin, zmax, m200min_SM, m200max_SM, size, seed=None):
         maximum value of the mass distribution, in units of solar masses
     size: int
         size of the distribution
+    seed: int
+        seed for random number generation
 
     Returns:
     -------

diff --git a/deepszsim/filters.py b/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
diff --git a/deepszsim/load_vars.py b/deepszsim/load_vars.py
@@ -10,7 +10,7 @@
 import sys
 import h5py
 
-def load_vars(file_name= "inputdata.yaml",
+def load_vars(file_name= "config_simACTDR5.yaml",
               file_dir = os.path.join(os.path.dirname(__file__), "Settings"),
               survey_num : int = None,
               survey_name : str = None,

diff --git a/deepszsim/make_sz_cluster.py b/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,10 @@ 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`
+        such as `500c`
 
     Returns:
     -------
@@ -439,9 +435,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
diff --git a/deepszsim/noise.py b/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.

diff --git a/deepszsim/read_yaml.py b/deepszsim/read_yaml.py
@@ -8,7 +8,7 @@
 
 class YAMLOperator:
 
-    def __init__(self, file_path=os.path.join(os.path.dirname(__file__), "Settings", "inputdata.yaml")):
+    def __init__(self, file_path=os.path.join(os.path.dirname(__file__), "Settings", "config_simACTDR5.yaml")):
 
         self.file_path = file_path