Synchrotron, new small scales

CHANGES.rst

@@ -1,7 +1,8 @@
 3.4.0 (unreleased)
 ==================
 
-- Brand new implementation of small scales injection for GNILC Dust `PR 133 <https://github.com/galsci/pysm/pull/133>`_, affects `d9, `d10`, `d11`
+- Brand new implementation of small scales injection for Synchrotron `PR 134 <https://github.com/galsci/pysm/pull/134>`_, affects `s4`, `s5`, `s6`, `s7`
+- Brand new implementation of small scales injection for GNILC Dust `PR 133 <https://github.com/galsci/pysm/pull/133>`_, affects `d9`, `d10`, `d11`
 - Fix bug in `InterpolatingComponent`, when the user requested a frequency between 2 available points, the weighting of the 2 relevant maps was switched, see `PR 129 <https://github.com/galsci/pysm/pull/129>`_
 - Implemented a proper unit test of the running `trapz` implementation used for bandpass integration against `np.trapz`, see `PR 129 <https://github.com/galsci/pysm/pull/129>`_
 - Imported WebSky extralactic components from `so_pysm_models`, now version 0.4, it also includes SPT based correction for CIB `PR 129 <https://github.com/galsci/pysm/pull/129>`_

docs/models.rst

@@ -29,8 +29,8 @@ The PySM 2 paper describes how input data (e.g. component separation maps from P
 most templates have been smoothed to remove high $\ell$ noise and then added small scale fluctuations.
 Here we try to reproduce the process to clarify it:
 
-* `Dust polarization templates from Planck Commander component separation outputs <preprocess-templates/reproduce_pysm2_dust_pol.html>`_, used in all dust models from `d0` to `d8` except `d4`
-* `Synchrotron polarization templates from WMAP low frequency maps <preprocess-templates/reproduce_pysm2_sync_pol.html>`_, used in all synchrotron models
+* `Dust polarization templates from Planck Commander component separation outputs <preprocess-templates/reproduce_pysm2_dust_pol.ipynb>`_, used in all dust models from `d0` to `d8` except `d4`
+* `Synchrotron polarization templates from WMAP low frequency maps <preprocess-templates/reproduce_pysm2_sync_pol.ipynb>`_, used in all synchrotron models
 
 
 Dust
@@ -54,7 +54,7 @@ Dust
 
 - **d9**: simplified version of **d10** with a fixed spectral index of 1.48 and a fixed dust black body temperature of 19.6 K all over the sky, based on Planck 2018 results.
 
-- **d10**: single component modified black body model based on templates from the `GNILC needlet-based analysis of Planck data <https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/Foreground_maps#GNILC_thermal_dust_maps>`_, with reduced contamination from CIB and point sources compared to the Commander maps used on **d1**. Small-scale fluctuations up to a $\ell_{max}$ of 16384 have been added to the templates in the `logpoltens` (Logarithm of the Polarization Fraction Tensor) formalism. Also the spectral index and dust temperature maps have random fluctuations at small scales. Available up to $N_{side}$ of 8192. Input templates are available `at NERSC <https://portal.nersc.gov/project/cmb/pysm-data/dust_gnilc/>`_.  Input templates are available at 2048, 4096 and 8192; lower nside are simulated by running `hp.ud_grade` on the $N_{side}=2048$ maps.  For more details about the pre-processing of the input data, see the `notebook about creating the GNILC-based template and adding small scales using the logpoltens formalism <preprocess-templates/small_scale_dust_pysm3.ipynb>`_ and `the notebook about spectral index and dust temperature <preprocess-templates/gnilc_dust_spectralindex_Tdust.html>`_. We also `apply a color correction <preprocess-templates/small_scale_dust_pysm3_generate_templates.ipynb>`_ to the 353 GHz template using the value provided by the `Planck 2018 XI paper <https://www.aanda.org/articles/aa/pdf/2020/09/aa32618-18.pdf>`_ (Table 2). The galactic region of about 3% of the sky within the GAL 97 Planck mask is set to the input GNILC map at 21.8' to avoid excess power in the spectra caused by the injected small scale signal.
+- **d10**: single component modified black body model based on templates from the `GNILC needlet-based analysis of Planck data <https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/Foreground_maps#GNILC_thermal_dust_maps>`_, with reduced contamination from CIB and point sources compared to the Commander maps used on **d1**. Small-scale fluctuations up to a $\ell_{max}$ of 16384 have been added to the templates in the `logpoltens` (Logarithm of the Polarization Fraction Tensor) formalism. Also the spectral index and dust temperature maps have random fluctuations at small scales. Available up to $N_{side}$ of 8192. Input templates are available `at NERSC <https://portal.nersc.gov/project/cmb/pysm-data/dust_gnilc/>`_.  Input templates are available at 2048, 4096 and 8192; lower nside are simulated by running `hp.ud_grade` on the $N_{side}=2048$ maps.  For more details about the pre-processing of the input data, see the `notebook about creating the GNILC-based template and adding small scales using the logpoltens formalism <preprocess-templates/small_scale_dust_pysm3.ipynb>`_ and `the notebook about spectral index and dust temperature <preprocess-templates/gnilc_dust_spectralindex_Tdust.ipynb>`_. We also `apply a color correction <preprocess-templates/small_scale_dust_pysm3_generate_templates.ipynb>`_ to the 353 GHz template using the value provided by the `Planck 2018 XI paper <https://www.aanda.org/articles/aa/pdf/2020/09/aa32618-18.pdf>`_ (Table 2). The galactic region of about 3% of the sky within the GAL 97 Planck mask is set to the input GNILC map at 21.8' to avoid excess power in the spectra caused by the injected small scale signal.
 
 - **d11**: like **d10** with stochastic small scales generated on-the-fly. It can reproduce **d10** if configured to run with a specific set of seeds and a specific $\ell_{max}$, see :py:class:`~pysm3.ModifiedBlackBodyRealization`. However reproducing **d10** is expensive because it needs to generate small scales with a $\ell_{max}=16384$ whatever output resolution is required, normally instead **d11** generates small scales just up to $\ell_{max}=2.5N_{side}$.
 
@@ -71,11 +71,11 @@ Synchrotron
 
 - **s4**: simplified version of **s5** with a fixed spectral index of -3.1 all over the sky.
 
-- **s5**: power law model based on the same templates of **s1**, Haslam in temperature and WMAP 9 year 23 GHz in polarization. Small-scale fluctuations up to a $\ell_{max}$ of 16384 have been added to the templates in the `logpoltens` (Logarithm of the Polarization Fraction Tensor) formalism. The spectral index map from **s1** has been rescaled `based on S-PASS <https://arxiv.org/abs/1802.01145>`_ and had small scales added to upgrade it up to $N_{side}$ 8192. Input templates are available `at NERSC <https://portal.nersc.gov/project/cmb/pysm-data/synch/>`_ at 2048, 4096 and 8192; lower nside are simulated by running `hp.ud_grade` on the $N_{side}=2048$ maps.  For more details about the pre-processing of the input data, see the `notebook about creating the templates and adding small scales using the logpoltens formalism <preprocess-templates/synchrotron_template_logpoltens.html>`_ and `the notebook about spectral index <preprocess-templates/synchrotron_beta.html>`_.
+- **s5**: power law model based on the same templates of **s1**, Haslam in temperature and WMAP 9 year 23 GHz in polarization. Small-scale fluctuations up to a $\ell_{max}$ of 16384 have been added to the templates in the `logpoltens` (Logarithm of the Polarization Fraction Tensor) formalism. The spectral index map from **s1** has been rescaled `based on S-PASS <https://arxiv.org/abs/1802.01145>`_ and had small scales added to upgrade it up to $N_{side}$ 8192. Input templates are available `at NERSC <https://portal.nersc.gov/project/cmb/pysm-data/synch/>`_ at 2048, 4096 and 8192; lower nside are simulated by running `hp.ud_grade` on the $N_{side}=2048$ maps.  For more details about the pre-processing of the input data, see the `notebook about creating the templates and adding small scales using the logpoltens formalism <preprocess-templates/synchrotron_template_logpoltens.ipynb>`_ and `the notebook about spectral index <preprocess-templates/synchrotron_beta.ipynb>`_.
 
 - **s6**: like **s5** with stochastic small scales generated on-the-fly. It can reproduce **s5** if configured to run with a specific set of seeds and a specific $\ell_{max}$, see :py:class:`~pysm3.PowerLawRealization`. However reproducing **s5** is expensive because it needs to generate small scales with a $\ell_{max}=16384$ whatever output resolution is required, normally instead **s6** generates small scales just up to $\ell_{max}=3N_{side}-1$.
 
-- **s7**: a power law with a curved index. The model uses the same templates and the same spectral index map of **s5**, the curvature term is based on the smoothed intensity template matched to the patch measured by the ARCADE experiment (Kogut, A. 2012, ApJ, 753, 110) and has random small scale fluctuations added, see `the relevant notebook <preprocess-templates/synchrotron_curvature.html>`_. The curvature map is available at $N_{side}=2048/4096/8192$.
+- **s7**: a power law with a curved index. The model uses the same templates and the same spectral index map of **s5**, the curvature term is based on the smoothed intensity template matched to the patch measured by the ARCADE experiment (Kogut, A. 2012, ApJ, 753, 110) and has random small scale fluctuations added, see `the relevant notebook <preprocess-templates/synchrotron_curvature.ipynb>`_. The curvature map is available at $N_{side}=2048/4096/8192$.
 
 
 AME

docs/preprocess-templates/utils_gnilc_generate_map_spectralindex_Td.ipynb

@@ -32,7 +32,11 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "tags": [
+     "parameters"
+    ]
+   },
    "outputs": [],
    "source": [
     "output_nside = 2048\n",
@@ -251,14 +255,22 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "pycmb",
+   "language": "python",
+   "name": "pycmb"
+  },
   "language_info": {
    "codemirror_mode": {
-    "name": "ipython"
+    "name": "ipython",
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
-   "nbconvert_exporter": "python"
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.0"
   }
  },
  "nbformat": 4,

docs/preprocess-templates/utils_synch_generate_map.ipynb

@@ -32,7 +32,7 @@
    },
    "outputs": [],
    "source": [
-    "output_nside = 8192"
+    "output_nside = 2048"
    ]
   },
   {
@@ -41,7 +41,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "output_lmax = min(3*output_nside - 1, 8192*2)"
+    "output_lmax = int(min(2.5*output_nside, 8192*2))"
    ]
   },
   {
@@ -67,7 +67,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "largescale_filename = list(output_dir_raw.glob(f\"{comp}*largescale*{sub}*.fits\"))"
+    "largescale_filename = list(output_dir_raw.glob(f\"{comp}*largescale*{sub}*.fits*\"))"
    ]
   },
   {
@@ -128,7 +128,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "modulate_alm = { k:hp.read_alm(output_dir_raw/f\"synch_{k}_modulation_alms_lmax1535.fits\").astype(np.complex128) for k in [\"temperature\",\"polarization\"] }"
+    "modulate_alm = { k:hp.read_alm(output_dir_raw/f\"synch_{k}_modulation_alms_lmax768.fits.gz\").astype(np.complex128) for k in [\"temperature\",\"polarization\"] }"
    ]
   },
   {
@@ -144,7 +144,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cl_small_scale = hp.read_cl(output_dir_raw / \"synch_small_scales_logpoltens_cl_lmax16384.fits\")"
+    "cl_small_scale = hp.read_cl(output_dir_raw / \"synch_small_scales_logpoltens_cl_lmax16384.fits.gz\")"
    ]
   },
   {
@@ -155,15 +155,16 @@
    "source": [
     "synalm_lmax = 8192*2 # for reproducibility\n",
     "# synalm_lmax = output_lmax\n",
+    "\n",
     "np.random.seed(555)\n",
     "\n",
     "alm_log_pol_tens_small_scale = hp.synalm(\n",
-    "    list(cl_small_scale) + [np.zeros_like(cl_small_scale[0])] * 3,\n",
+    "    list(cl_small_scale),\n",
     "    lmax=synalm_lmax,\n",
     "    new=True,\n",
     ")\n",
     "\n",
-    "alm_log_pol_tens_small_scale = [hp.almxfl(each, np.ones(3*output_nside-1)) for each in alm_log_pol_tens_small_scale]\n",
+    "alm_log_pol_tens_small_scale = [hp.almxfl(each, np.ones(output_lmax+1)) for each in alm_log_pol_tens_small_scale]\n",
     "map_log_pol_tens_small_scale = hp.alm2map(alm_log_pol_tens_small_scale, nside=output_nside)\n",
     "map_log_pol_tens_small_scale[0] *= hp.alm2map(modulate_alm[\"temperature\"], output_nside)\n",
     "map_log_pol_tens_small_scale[1:] *= hp.alm2map(modulate_alm[\"polarization\"], output_nside)\n",
@@ -225,20 +226,13 @@
    "source": [
     "add_metadata([output_dir / f\"synch_template_nside{output_nside}.fits\"], coord=\"G\", unit=\"uK_RJ\", ref_freq=\"23 GHz\")"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "condanamaster2",
+   "display_name": "pycmb",
    "language": "python",
-   "name": "condanamaster2"
+   "name": "pycmb"
   },
   "language_info": {
    "codemirror_mode": {
@@ -250,7 +244,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.0"
+   "version": "3.10.0"
   }
  },
  "nbformat": 4,

docs/preprocess-templates/utils_synch_generate_map_spectralindex.ipynb

@@ -41,7 +41,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "output_lmax = min(3 * output_nside - 1, 8192 * 2)"
+    "output_lmax = int(min(2.5 * output_nside, 8192 * 2))"
    ]
   },
   {
@@ -58,7 +58,7 @@
    "outputs": [],
    "source": [
     "alm_large_scale = hp.read_alm(\n",
-    "    output_dir_raw / f\"synch_largescale_beta_alm_nside512_lmax1535_complex64.fits\",\n",
+    "    output_dir_raw / f\"synch_largescale_beta_alm_nside512_lmax1535_complex64.fits.gz\",\n",
     "    hdu=1,\n",
     ")"
    ]
@@ -86,7 +86,7 @@
    "outputs": [],
    "source": [
     "modulate_alm = hp.read_alm(\n",
-    "    output_dir_raw / f\"synch_temperature_modulation_alms_lmax1535.fits\"\n",
+    "    output_dir_raw / f\"synch_temperature_modulation_alms_lmax768.fits.gz\"\n",
     ").astype(np.complex128)"
    ]
   },
@@ -104,7 +104,7 @@
    "outputs": [],
    "source": [
     "cl_small_scale = hp.read_cl(\n",
-    "    output_dir_raw / f\"synch_small_scales_beta_cl_lmax16384.fits\"\n",
+    "    output_dir_raw / f\"synch_small_scales_beta_cl_lmax16384.fits.gz\"\n",
     ")"
    ]
   },
@@ -125,7 +125,7 @@
     "    new=True,\n",
     ")\n",
     "\n",
-    "alm_small_scale = hp.almxfl(alm_small_scale, np.ones(3 * output_nside - 1))\n",
+    "alm_small_scale = hp.almxfl(alm_small_scale, np.ones(output_lmax+1))\n",
     "map_small_scale = hp.alm2map(alm_small_scale, nside=output_nside)\n",
     "map_small_scale *= hp.alm2map(modulate_alm, output_nside)\n",
     "assert np.isnan(map_small_scale).sum() == 0"
@@ -198,9 +198,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "condanamaster2",
+   "display_name": "pycmb",
    "language": "python",
-   "name": "condanamaster2"
+   "name": "pycmb"
   },
   "language_info": {
    "codemirror_mode": {
@@ -212,7 +212,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.0"
+   "version": "3.10.0"
   }
  },
  "nbformat": 4,

pysm3/data/presets.cfg

@@ -245,12 +245,12 @@ freq_ref_P = "23 GHz"
 max_nside = 8192
 [s6]
 class = "PowerLawRealization"
-largescale_alm = "synch/raw/synch_largescale_template_logpoltens_alm_nside512_lmax1535_complex64.fits"
-amplitude_modulation_temp_alm = "synch/raw/synch_temperature_modulation_alms_lmax1535.fits"
-amplitude_modulation_pol_alm = "synch/raw/synch_polarization_modulation_alms_lmax1535.fits"
-small_scale_cl = "synch/raw/synch_small_scales_logpoltens_cl_lmax16384.fits"
-largescale_alm_pl_index = "synch/raw/synch_largescale_beta_alm_nside512_lmax1535_complex64.fits"
-small_scale_cl_pl_index = "synch/raw/synch_small_scales_beta_cl_lmax16384.fits"
+largescale_alm = "synch/raw/synch_largescale_template_logpoltens_alm_nside512_lmax768_complex64.fits.gz"
+amplitude_modulation_temp_alm = "synch/raw/synch_temperature_modulation_alms_lmax768.fits.gz"
+amplitude_modulation_pol_alm = "synch/raw/synch_polarization_modulation_alms_lmax768.fits.gz"
+small_scale_cl = "synch/raw/synch_small_scales_logpoltens_cl_lmax16384.fits.gz"
+largescale_alm_pl_index = "synch/raw/synch_largescale_beta_alm_nside512_lmax1535_complex64.fits.gz"
+small_scale_cl_pl_index = "synch/raw/synch_small_scales_beta_cl_lmax16384.fits.gz"
 freq_ref = "23 GHz"
 max_nside = 8192
 # Configuration for reproducing s5

pysm3/models/power_law_realization.py

@@ -96,19 +96,20 @@ def draw_realization(self, synalm_lmax=None, seeds=None):
             seeds = (None, None)
 
         if synalm_lmax is None:
-            synalm_lmax = min(16384, 3 * self.nside - 1)
+            synalm_lmax = int(min(16384, 2.5 * self.nside))
+
+        output_lmax = int(min(synalm_lmax, 2.5 * self.nside))
 
         np.random.seed(seeds[0])
 
         alm_small_scale = hp.synalm(
-            list(self.small_scale_cl.value)
-            + [np.zeros_like(self.small_scale_cl[0])] * 3,
+            list(self.small_scale_cl.value),
             lmax=synalm_lmax,
             new=True,
         )
 
         alm_small_scale = [
-            hp.almxfl(each, np.ones(min(synalm_lmax, 3 * self.nside - 1)))
+            hp.almxfl(each, np.ones(output_lmax+1))
             for each in alm_small_scale
         ]
         map_small_scale = hp.alm2map(alm_small_scale, nside=self.nside)
@@ -138,7 +139,7 @@ def draw_realization(self, synalm_lmax=None, seeds=None):
         )
 
         alm_small_scale = hp.almxfl(
-            alm_small_scale, np.ones(min(3 * self.nside - 1, synalm_lmax + 1))
+            alm_small_scale, np.ones(output_lmax + 1)
         )
         pl_index = hp.alm2map(alm_small_scale, nside=self.nside) * output_unit
         pl_index *= modulate_map_I

pysm3/tests/test_synch_pysm3.py

@@ -115,10 +115,11 @@ def test_synch_44(model_tag):
     psutil.virtual_memory().total * u.byte < 20 * u.GB,
     reason="Running s6 at high lmax requires 20 GB of RAM",
 )
-def test_s6_vs_s5():
+@pytest.mark.parametrize("freq", [23, 44])
+def test_s6_vs_s5(freq):
     nside = 2048
 
-    freq = 44 * u.GHz
+    freq = freq * u.GHz
 
     output_s5 = pysm3.Sky(preset_strings=["s5"], nside=nside).get_emission(freq)
     s6_configuration = pysm3.sky.PRESET_MODELS["s6"].copy()