diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml index 84b691e4..08198d19 100644 --- a/.github/workflows/ci_tests.yml +++ b/.github/workflows/ci_tests.yml @@ -37,18 +37,43 @@ jobs: uses: actions/checkout@v4 with: fetch-depth: 0 - - uses: mamba-org/setup-micromamba@v1 + - name: Free disk space + run: | + sudo rm -rf /usr/share/dotnet /usr/local/lib/android /opt/ghc /opt/hostedtoolcache + df -h + shell: bash -el {0} + - name: Prepare large temp dir + run: | + if [ -d /mnt ]; then + sudo mkdir -p /mnt/runner-temp + sudo chown "$USER":"$USER" /mnt/runner-temp + echo "PYSM_TEMP_ROOT=/mnt/runner-temp" >> "$GITHUB_ENV" + else + echo "PYSM_TEMP_ROOT=$RUNNER_TEMP" >> "$GITHUB_ENV" + fi + shell: bash -el {0} + - name: Set temp locations + run: | + echo "ASTROPY_CACHE_DIR=$PYSM_TEMP_ROOT/astropy-cache" >> "$GITHUB_ENV" + echo "ASTROPY_CONFIG_DIR=$PYSM_TEMP_ROOT/astropy-config" >> "$GITHUB_ENV" + echo "XDG_CACHE_HOME=$PYSM_TEMP_ROOT/xdg-cache" >> "$GITHUB_ENV" + echo "XDG_CONFIG_HOME=$PYSM_TEMP_ROOT/xdg-config" >> "$GITHUB_ENV" + echo "TMPDIR=$PYSM_TEMP_ROOT/astropy-tmp" >> "$GITHUB_ENV" + echo "HOME=$PYSM_TEMP_ROOT/home" >> "$GITHUB_ENV" + shell: bash -el {0} + - name: Prepare cache directories + run: | + mkdir -p "$HOME" "$ASTROPY_CACHE_DIR" "$ASTROPY_CONFIG_DIR" \ + "$XDG_CACHE_HOME" "$XDG_CONFIG_HOME" "$TMPDIR" + shell: bash -el {0} + - name: Set up Python + uses: actions/setup-python@v5 with: - environment-name: test-env - create-args: >- - python=${{ matrix.python }} - netcdf4 - init-shell: >- - bash - cache-environment: true - post-cleanup: 'all' + python-version: ${{ matrix.python }} - name: Install package - run: pip install .[test] + run: | + python -m pip install --upgrade pip + python -m pip install --no-cache-dir .[test] shell: bash -el {0} - name: Run Pytest run: | diff --git a/.github/workflows/import_base_requirements.yaml b/.github/workflows/import_base_requirements.yaml index 3b810264..ba49a22b 100644 --- a/.github/workflows/import_base_requirements.yaml +++ b/.github/workflows/import_base_requirements.yaml @@ -22,20 +22,40 @@ jobs: uses: actions/checkout@v4 with: fetch-depth: 0 - - uses: mamba-org/setup-micromamba@v1 + - name: Prepare large temp dir + run: | + if [ -d /mnt ]; then + sudo mkdir -p /mnt/runner-temp + sudo chown "$USER":"$USER" /mnt/runner-temp + echo "PYSM_TEMP_ROOT=/mnt/runner-temp" >> "$GITHUB_ENV" + else + echo "PYSM_TEMP_ROOT=$RUNNER_TEMP" >> "$GITHUB_ENV" + fi + shell: bash -el {0} + - name: Set temp locations + run: | + echo "ASTROPY_CACHE_DIR=$PYSM_TEMP_ROOT/astropy-cache" >> "$GITHUB_ENV" + echo "ASTROPY_CONFIG_DIR=$PYSM_TEMP_ROOT/astropy-config" >> "$GITHUB_ENV" + echo "XDG_CACHE_HOME=$PYSM_TEMP_ROOT/xdg-cache" >> "$GITHUB_ENV" + echo "XDG_CONFIG_HOME=$PYSM_TEMP_ROOT/xdg-config" >> "$GITHUB_ENV" + echo "TMPDIR=$PYSM_TEMP_ROOT/astropy-tmp" >> "$GITHUB_ENV" + echo "HOME=$PYSM_TEMP_ROOT/home" >> "$GITHUB_ENV" + shell: bash -el {0} + - name: Prepare cache directories + run: | + mkdir -p "$HOME" "$ASTROPY_CACHE_DIR" "$ASTROPY_CONFIG_DIR" \ + "$XDG_CACHE_HOME" "$XDG_CONFIG_HOME" "$TMPDIR" + shell: bash -el {0} + - name: Set up Python + uses: actions/setup-python@v5 with: - environment-name: test-env - create-args: >- - python=${{ matrix.python }} - netcdf4 - init-shell: >- - bash - cache-environment: true - post-cleanup: 'all' + python-version: ${{ matrix.python }} - name: Install base requirements - run: pip install . + run: | + python -m pip install --upgrade pip + python -m pip install --no-cache-dir . shell: bash -el {0} - name: Check import run: | python -c "import pysm3" - shell: bash -el {0} \ No newline at end of file + shell: bash -el {0} diff --git a/CHANGES.rst b/CHANGES.rst index fd773fef..e4cc6ebe 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,7 @@ Unreleased ========== - Add NumPy 2 support by replacing deprecated `np.trapz` with `trapezoid` https://github.com/galsci/pysm/pull/232 +- Add HalfDome TSZ maps (new 'tsz2' model) to pysm. Create new documentation and code following the same structural approach used for WebSky. 3.4.3 (2025-10-02) ================== diff --git a/docs/halfdome.rst b/docs/halfdome.rst new file mode 100644 index 00000000..7379fa3a --- /dev/null +++ b/docs/halfdome.rst @@ -0,0 +1,37 @@ +.. _halfdome: + +HalfDome +****** + +HalfDome cosmological simulations are a set of simulations tailored for the joint analysis of Stage-IV surveys. For more information, see: [https://halfdomesims.github.io/] (https://halfdomesims.github.io/) + +Cosmic Infrared Background +-------------------------- + +TODO + +Radio Galaxies +-------------- + +TODO + +Thermal SZ Effect +----------------- + +Provided are eleven realizations of the Compton-y parameter. They were generated with xgpaint bataglia16 profiles at $N_{side}=8192$. + +Kinetic SZ Effect +----------------- + +TODO + + +Lensing Convergence +------------------- + +TODO + +Primary and lensed CMB +---------------------- + +TODO diff --git a/docs/models.rst b/docs/models.rst index 1fc12c3d..e2e780db 100644 --- a/docs/models.rst +++ b/docs/models.rst @@ -138,8 +138,11 @@ Sunyaev–Zeldovich emission - **tsz1**: Thermal SZ emission from WebSky 0.4. Available at $N_{side}=8192$. For more details see :ref:`websky`. +- **tsz2**: Thermal SZ emission from HalfDome 0.1. Generated using xgpaint with Battaglia16 profiles. Eleven realizations available at $N_{side}=8192$. For more details see :ref:`halfdome`. + - **ksz1**: Kinetic SZ emission from WebSky 0.4. Available at $N_{side}=4096$. For more details see :ref:`websky`. + Radio galaxies ============== diff --git a/src/pysm3/data/presets.cfg b/src/pysm3/data/presets.cfg index 006b1ece..ec82881f 100644 --- a/src/pysm3/data/presets.cfg +++ b/src/pysm3/data/presets.cfg @@ -394,6 +394,13 @@ class = "WebSkySZ" version = "0.4" sz_type = "thermal" max_nside = 4096 +[tsz2] +class = "HalfDomeSZ" +version = "0.1" +sz_type = "thermal" +map_sz = "y_b16_halo_res1_s{seed}.fits" +seed = 0 +max_nside = 8192 [dip1] # From NPIPE https://doi.org/10.1051/0004-6361/202038073 # diff --git a/src/pysm3/models/__init__.py b/src/pysm3/models/__init__.py index 44dbbb56..a1529d77 100644 --- a/src/pysm3/models/__init__.py +++ b/src/pysm3/models/__init__.py @@ -18,5 +18,13 @@ WebSkySZ, WebSkyRadioGalaxies, ) +from .halfdome import ( + HalfDomeSZ, +# SPT_CIB_map_scaling, +# HalfDomeCMB, +# HalfDomeCIB, +# HalfDomeSZ, +# HalfDomeRadioGalaxies, +) from .catalog import PointSourceCatalog from .dipole import CMBDipole diff --git a/src/pysm3/models/catalog.py b/src/pysm3/models/catalog.py index 5c456385..ea96ca59 100644 --- a/src/pysm3/models/catalog.py +++ b/src/pysm3/models/catalog.py @@ -152,11 +152,16 @@ def __init__( else: self.catalog_slice = catalog_slice + def _attr_to_str(value): + if isinstance(value, (bytes, np.bytes_)): + return value.decode("utf-8") + return str(value) + with h5py.File(self.catalog_filename) as f: - assert f["theta"].attrs["units"].decode("UTF-8") == "rad" - assert f["phi"].attrs["units"].decode("UTF-8") == "rad" - assert f["logpolycoefflux"].attrs["units"].decode("UTF-8") == "Jy" - assert f["logpolycoefpolflux"].attrs["units"].decode("UTF-8") == "Jy" + assert _attr_to_str(f["theta"].attrs["units"]) == "rad" + assert _attr_to_str(f["phi"].attrs["units"]) == "rad" + assert _attr_to_str(f["logpolycoefflux"].attrs["units"]) == "Jy" + assert _attr_to_str(f["logpolycoefpolflux"].attrs["units"]) == "Jy" assert map_dist is None, "Distributed execution not supported" diff --git a/src/pysm3/models/halfdome.py b/src/pysm3/models/halfdome.py new file mode 100644 index 00000000..dbf521fd --- /dev/null +++ b/src/pysm3/models/halfdome.py @@ -0,0 +1,323 @@ +import os.path +from pathlib import Path + +import numpy as np +from numba import njit + +import pysm3 as pysm +import pysm3.units as u + +from .. import utils +from .cmb import CMBMap +from .interpolating import InterpolatingComponent +from .template import Model + +# Seeds for the HalfDome simulations. +SEEDS = [100, 102, 104, 106, 108, 110, 112, 114, 116, 118, 120] + + +@njit +def y2uK_CMB(nu): + """Compton-y distortion at a given frequency + Parmeters: + nu (float): frequency in GHz + Returns: + float: intensity variation dT_CMB in micro-Kelvin + dT_CMB = dI_nu / (dB_nu / dT)_Tcmb + where B_nu is the Planck function and dI_nu is the intensity distortion + """ + + h = 6.62607004e-27 + k = 1.380622e-16 + Tcmb = 2.725 + x = h * nu * 1e9 / k / Tcmb + return 1e6 * Tcmb * (x * (np.exp(x) + 1) / (np.exp(x) - 1) - 4) + + +class HalfDomeSZ(Model): + def __init__( + self, + nside, + version="0.1", + sz_type="thermal", + map_sz="y_b16_halo_res1_s{seed}.fits", + seed = 0, + max_nside=None, + map_dist=None, + ): + + if max_nside is None: + if sz_type == "kinetic": + # max_nside = 4096 + raise NotImplementedError("Kinetic SZ maps are not implemented yet.") + if sz_type == "thermal": + max_nside = 8192 + super().__init__(nside=nside, max_nside=max_nside, map_dist=map_dist) + self.version = str(version) + self.sz_type = sz_type + self.map_sz = map_sz.format(seed=SEEDS[seed]) + self.remote_data = utils.RemoteData() + filename = self.remote_data.get(self.get_filename()) + self.m = self.read_map(filename, field=0, unit=u.uK_CMB) + + def get_filename(self): + """Get SZ filenames for a halfdome version""" + + path = Path("halfdome") / self.version + if self.sz_type == "kinetic": + path = path / "ksz" / self.map_sz + elif self.sz_type == "thermal": + path = path / "tsz" / self.map_sz + + return str(path) + + @u.quantity_input + def get_emission( + self, freqs: u.Quantity[u.GHz], weights=None + ) -> u.Quantity[u.uK_RJ]: + + freqs = pysm.check_freq_input(freqs) + weights = pysm.normalize_weights(freqs, weights) + + # input map is in uK_CMB, we multiply the weights which are + # in uK_RJ by the conversion factor of uK_CMB->uK_RJ + # this is the equivalent of + weights = (weights * u.uK_CMB).to_value( + u.uK_RJ, equivalencies=u.cmb_equivalencies(freqs * u.GHz) + ) + + is_thermal = self.sz_type == "thermal" + return ( + get_sz_emission_numba(freqs, weights, self.m.value, is_thermal) << u.uK_RJ + ) + + # the output of out is always 2D, (IQU, npix) + + +@njit(parallel=True) +def get_sz_emission_numba(freqs, weights, m, is_thermal): + output = np.zeros((3, len(m)), dtype=np.float64) + for i in range(len(freqs)): + signal = m * y2uK_CMB(freqs[i]) if is_thermal else m.astype(np.float64) + pysm.utils.trapz_step_inplace(freqs, weights, i, signal, output[0]) + return output + +''' +@njit +def SPT_CIB_map_scaling(nu): + """CIB maps correction based on SPT data + + Parameters + ---------- + nu : float or np.array + frequency in GHz + + Returns + ------- + correction : float or np.array + correction factor to be applied to the map + """ + return np.sqrt(1 + (1.84 - 1) / (1 + np.exp((nu - 227) / 75))) + + +class HalfDomeCIB(InterpolatingComponent): + """TODO PySM component interpolating between precomputed maps""" + + def __init__( + self, + nside, + halfdome_version="0.1", + input_units="MJy / sr", + max_nside=4096, + interpolation_kind="linear", + apply_SPT_correction=True, + local_folder=None, + map_dist=None, + ): + """Load and interpolate HalfDome CIB maps + + Parameters + ---------- + nside : nside + target nside of the output maps + halfdome_version : str + currently only 0.1 is supported + input_units : str + input units string, e.g. uK_CMB, K_RJ + max_nside : int + maximum nside at which the input maps are available at + `nside` can be higher than this, but then PySM will use + `ud_grade` to create maps at higher resolution. + interpolation_kind : str + See the docstring of :py:class:`~pysm3.InterpolatingComponent` + apply_SPT_correction : bool + Apply the correction computed by comparison with the South + Pole Telescope maps. + local_folder : str + Override the input maps folder + map_dist : :py:class:`~pysm3.MapDistribution` + Required for partial sky or MPI, see the PySM docs + """ + self.local_folder = local_folder + self.halfdome_freqs_float = [ + 17.0, + 18.7, + 21.6, + 24.5, + 27.3, + 30.0, + 35.9, + 41.7, + 44.0, + 47.4, + 63.9, + 67.8, + 70.0, + 73.7, + 79.6, + 90.2, + 100, + 111, + 129, + 143, + 153, + 164, + 189, + 210, + 217, + 232, + 256, + 275, + 294, + 306, + 314, + 340, + 353, + 375, + 409, + 467, + 525, + 545, + 584, + 643, + 729, + 817, + 857, + 906, + 994, + 1080, + ] + self.halfdome_freqs = [f"{f:06.1f}" for f in self.halfdome_freqs_float] + self.apply_SPT_correction = apply_SPT_correction + super().__init__( + path=halfdome_version, + input_units=input_units, + nside=nside, + max_nside=max_nside, + interpolation_kind=interpolation_kind, + map_dist=map_dist, + ) + self.remote_data = utils.RemoteData() + + def get_filenames(self, path): + """Get filenames for a halfdome version + For a standard interpolating component, we list files in folder, + here we need to know the names in advance so that we can only download + the required maps. + """ + + halfdome_version = path + + filenames = { + float(str_freq): f"halfdome/{halfdome_version}/cib/cib_{str_freq}.fits" + for str_freq in self.halfdome_freqs + } + + if self.local_folder is not None: + for freq in filenames: + filenames[freq] = os.path.join(self.local_folder, filenames[freq]) + return filenames + + def read_map_by_frequency(self, freq): + filename = self.remote_data.get(self.maps[freq]) + m = self.read_map_file(freq, filename) + if self.apply_SPT_correction: + m *= SPT_CIB_map_scaling(freq) + return m + + +# radio galaxies are just like CIB, just interpolating +class HalfDomeRadioGalaxies(HalfDomeCIB): + def __init__( + self, + nside, + halfdome_version="0.1", + input_units="Jy / sr", + max_nside=4096, + interpolation_kind="linear", + apply_SPT_correction=False, + local_folder=None, + map_dist=None, + ): + super().__init__( + nside=nside, + halfdome_version=halfdome_version, + input_units=input_units, + max_nside=max_nside, + interpolation_kind=interpolation_kind, + apply_SPT_correction=apply_SPT_correction, + local_folder=local_folder, + map_dist=map_dist, + ) + + def get_filenames(self, path): + """Get filenames for a halfdome version + For a standard interpolating component, we list files in folder, + here we need to know the names in advance so that we can only download + the required maps. + """ + + halfdome_version = path + + filenames = { + float(str_freq): f"halfdome/{halfdome_version}/radio/radio_{str_freq}.fits" + for str_freq in self.halfdome_freqs + } + + if self.local_folder is not None: + for freq in filenames: + filenames[freq] = os.path.join(self.local_folder, filenames[freq]) + return filenames + + +class HalfDomeCMB(CMBMap): + def __init__( + self, + nside, + max_nside=4096, + halfdome_version=0.1, + seed=1, + lensed=True, + include_solar_dipole=False, + map_dist=None, + ): + template_nside = 512 if nside <= 512 else 4096 + lens = "" if lensed else "un" + soldip = "solardipole_" if include_solar_dipole else "" + filenames = [ + utils.RemoteData().get( + f"halfdome/{halfdome_version}/cmb/map_{pol}_{lens}" + + f"lensed_alm_seed{seed}_{soldip}nside{template_nside}.fits" + ) + for pol in "IQU" + ] + super().__init__( + map_I=filenames[0], + map_Q=filenames[1], + map_U=filenames[2], + nside=nside, + max_nside=max_nside, + map_dist=map_dist, + ) + +''' diff --git a/tests/test_catalog.py b/tests/test_catalog.py index 3da90608..537863e0 100644 --- a/tests/test_catalog.py +++ b/tests/test_catalog.py @@ -1,9 +1,7 @@ -# Important: it is important to import netCDF4 before `h4py` is imported -# to avoid "HDF Error" under Ubuntu with pip -# This does not happen with conda packages +# Important: import netCDF4 before h5py to avoid "HDF Error" under Ubuntu with pip. +import netCDF4 import h5py import healpy as hp -import netCDF4 import numpy as np try: diff --git a/tests/test_halfdome.py b/tests/test_halfdome.py new file mode 100644 index 00000000..7f87b0c9 --- /dev/null +++ b/tests/test_halfdome.py @@ -0,0 +1,75 @@ +import os + +import healpy as hp +import numpy as np +import pytest + +try: # PySM >= 3.2.1 + import pysm3.units as u +except ImportError: + import pysm.units as u + +from pysm3 import ( + HalfDomeSZ, + # SPT_CIB_map_scaling, + # HalfDomeCIB, + # HalfDomeRadioGalaxies, + utils, +) # , WebSkyCMBTensor + + +@pytest.mark.parametrize("sz_type", ["thermal"]) # , "kinetic"]) (kinetic not implemented yet) +def test_halfdome_sz_seeds(tmp_path, monkeypatch, sz_type): + + os.environ.pop("PYSM_LOCAL_DATA", None) + monkeypatch.setattr(utils.data, "PREDEFINED_DATA_FOLDERS", [str(tmp_path)]) + + nside = 8 + shape = hp.nside2npix(nside) + path = tmp_path / "halfdome" / "0.1" / "tsz" + path.mkdir(parents=True) + hp.write_map(path / "y_b16_halo_res1_s100.fits", np.ones(shape, dtype=np.float32)) + hp.write_map( + path / "y_b16_halo_res1_s102.fits", np.ones(shape, dtype=np.float32) * 2 + ) + + tsz_0 = HalfDomeSZ(nside, "0.1", sz_type=sz_type, seed=0) + tsz_1 = HalfDomeSZ(nside, "0.1", sz_type=sz_type, seed=1) + + freq = 100 * u.GHz + tsz_0_map = tsz_0.get_emission(freq).to( + u.uK_CMB, equivalencies=u.cmb_equivalencies(freq) + ) + tsz_1_map = tsz_1.get_emission(freq).to( + u.uK_CMB, equivalencies=u.cmb_equivalencies(freq) + ) + + with pytest.raises(AssertionError): + np.testing.assert_allclose(tsz_0_map, tsz_1_map) + + +@pytest.mark.parametrize("sz_type", ["thermal"]) # , "kinetic"]) (kinetic not implemented yet) +def test_halfdome_sz(tmp_path, monkeypatch, sz_type): # It needs to be the last + + os.environ.pop("PYSM_LOCAL_DATA", None) + monkeypatch.setattr(utils.data, "PREDEFINED_DATA_FOLDERS", [str(tmp_path)]) + nside = 4 + shape = hp.nside2npix(nside) + + path = tmp_path / "halfdome" / "0.1" / "tsz" + path.mkdir(parents=True) + filename = "y_b16_halo_res1_s100.fits" if sz_type == "thermal" else "ksz.fits" + test_map = np.ones(shape, dtype=np.float32) + if sz_type == "thermal": + test_map *= 1e-6 + hp.write_map(path / filename, test_map) + + tsz = HalfDomeSZ(nside, "0.1", sz_type=sz_type, seed=0) + + freq = 100 * u.GHz + tsz_map = tsz.get_emission(freq).to( + u.uK_CMB, equivalencies=u.cmb_equivalencies(freq) + ) + value = -4.109055 * u.uK_CMB if sz_type == "thermal" else 1.0 * u.uK_CMB + np.testing.assert_allclose(np.ones(len(tsz_map[0])) * value, tsz_map[0], rtol=1e-4) + np.testing.assert_allclose(np.zeros((2, len(tsz_map[0]))) * u.uK_CMB, tsz_map[1:])