Synchrotron with local-spectra modulation

docs/preprocess-templates/utils_synch_generate_map.ipynb

@@ -115,6 +115,15 @@
     "map_log_pol_tens_large_scale = hp.alm2map(alm_log_pol_tens_large_scale.astype(np.complex128), nside=output_nside)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "del alm_log_pol_tens_large_scale"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -128,7 +137,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "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\"] }"
+    "modulate_alm = { k:hp.read_alm(output_dir_raw/f\"synch_{k}_modulation_alms_lmax192.fits.gz\").astype(np.complex128) for k in [\"temperature\",\"polarization\"] }"
    ]
   },
   {
@@ -144,7 +153,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cl_small_scale = hp.read_cl(output_dir_raw / \"synch_small_scales_logpoltens_cl_lmax16384.fits.gz\")"
+    "cl_small_scale = hp.read_cl(output_dir_raw / \"synch_small_scales_cl_lmax16384.fits.gz\")"
    ]
   },
   {
@@ -154,7 +163,7 @@
    "outputs": [],
    "source": [
     "synalm_lmax = 8192*2 # for reproducibility\n",
-    "# synalm_lmax = output_lmax\n",
+    "# synalm_lmax = 512\n",
     "\n",
     "np.random.seed(555)\n",
     "\n",
@@ -166,6 +175,7 @@
     "\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",
+    "del alm_log_pol_tens_small_scale\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",
     "assert np.isnan(map_log_pol_tens_small_scale).sum() == 0"
@@ -191,6 +201,15 @@
     "map_log_pol_tens = map_log_pol_tens_large_scale + map_log_pol_tens_small_scale"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "del map_log_pol_tens_large_scale, map_log_pol_tens_small_scale"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -209,6 +228,50 @@
     "output_map = log_pol_tens_to_map(map_log_pol_tens)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Galactic plane fix"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "galplane_fix = hp.read_map(output_dir_raw / \"synch_galplane.fits.gz\", (0, 1, 2, 3))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "output_map *= hp.ud_grade(galplane_fix[3], output_nside)\n",
+    "output_map += hp.ud_grade(galplane_fix[:3] * (1 - galplane_fix[3]), output_nside)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if output_nside < 4096:\n",
+    "    hp.mollview((output_map - log_pol_tens_to_map(map_log_pol_tens))[0])\n",
+    "    hp.mollview(output_map[0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Write outputs"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -226,6 +289,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": {

docs/preprocess-templates/utils_synch_generate_map.slurm

@@ -0,0 +1,25 @@
+#!/bin/bash
+#SBATCH --qos=regular
+#SBATCH --time=6:00:00
+#SBATCH --nodes=1
+#SBATCH --tasks-per-node=1
+#SBATCH --constraint=knl
+#SBATCH --error=synch_generate_map_%j.err
+#SBATCH --output=synch_generate_map_%j.out
+#SBATCH [email protected]
+#SBATCH --mail-type=ALL
+
+export PYTHONUNBUFFERED=1
+export OMP_PROC_BIND=true
+export OMP_PLACES=threads
+export OMP_NUM_THREADS=68
+export HDF5_USE_FILE_LOCKING=FALSE
+
+PAPERMILL="/global/common/software/cmb/zonca/conda/pycmb/bin/papermill"
+for NSIDE in 4096 8192
+do
+    for NB in "utils_synch_generate_map.ipynb"
+    do
+        srun --cpu_bind=cores $PAPERMILL $NB data/${NB/.ipynb/.$NSIDE.$NAME.ipynb} -p output_nside $NSIDE -p name $NAME
+    done
+done

docs/preprocess-templates/utils_synch_generate_map_curvature.ipynb

@@ -37,7 +37,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "output_lmax = min(3 * output_nside - 1, 8192 * 2)"
+    "output_lmax = min(int(2.5 * output_nside), 8192 * 2)"
    ]
   },
   {
@@ -54,7 +54,7 @@
    "outputs": [],
    "source": [
     "alm = hp.read_alm(\n",
-    "    output_dir_raw / f\"synch_curvature_alm_nside8192_lmax16384_complex64.fits\",\n",
+    "    output_dir_raw / f\"synch_curvature_alm_nside8192_lmax16384.fits\",\n",
     "    hdu=1,\n",
     ")"
    ]
@@ -111,19 +111,27 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "add_metadata([output_filename], coord=\"G\", unit=\"\", ref_freq=\"23 GHz\")"
+    "add_metadata([output_filename], coord=\"G\", unit=\"\", ref_freq=\"23 GHz\", ell_max=str(output_lmax))"
    ]
   }
  ],
  "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_spectralindex.ipynb

@@ -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.gz\",\n",
+    "    output_dir_raw / f\"synch_largescale_beta_alm_nside512_lmax768.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_lmax768.fits.gz\"\n",
+    "    output_dir_raw / f\"synch_amplitude_modulation_alms_lmax768.fits.gz\"\n",
     ").astype(np.complex128)"
    ]
   },
@@ -174,7 +174,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "add_metadata([output_dir / f\"synch_beta_nside{output_nside}.fits\"], coord=\"G\", unit=\"\", ref_freq=\"23 GHz\")"
+    "add_metadata([output_dir / f\"synch_beta_nside{output_nside}.fits\"], coord=\"G\", unit=\"\", ref_freq=\"23 GHz\", ell_max=str(output_lmax))"
    ]
   },
   {
@@ -187,13 +187,6 @@
     "hp.mollview(map_small_scale, title=\"Small scale\")\n",
     "hp.mollview(output_map, title=\"Total\")"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

pysm3/data/presets.cfg

@@ -246,14 +246,17 @@ freq_ref_P = "23 GHz"
 max_nside = 8192
 [s6]
 class = "PowerLawRealization"
-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"
+largescale_alm = "synch/raw/synch_largescale_template_logpoltens_alm_nside512_lmax1536.fits.gz"
+amplitude_modulation_temp_alm = "synch/raw/synch_temperature_modulation_alms_lmax192.fits.gz"
+amplitude_modulation_pol_alm = "synch/raw/synch_polarization_modulation_alms_lmax192.fits.gz"
+amplitude_modulation_beta_alm = "synch/raw/synch_amplitude_modulation_alms_lmax768.fits.gz"
+small_scale_cl = "synch/raw/synch_small_scales_cl_lmax16384.fits.gz"
+largescale_alm_pl_index = "synch/raw/synch_largescale_beta_alm_nside512_lmax768.fits.gz"
 small_scale_cl_pl_index = "synch/raw/synch_small_scales_beta_cl_lmax16384.fits.gz"
 freq_ref = "23 GHz"
 max_nside = 8192
+# Remove the galplane_fix option to skip applying the galactic plane fix
+galplane_fix = "synch/raw/synch_galplane.fits.gz"
 # Configuration for reproducing s5
 # seeds = [555,444]
 # synalm_lmax = 16384

pysm3/models/power_law_realization.py

@@ -18,6 +18,8 @@ def __init__(
         largescale_alm_pl_index,
         small_scale_cl_pl_index,
         nside,
+        amplitude_modulation_beta_alm=None,
+        galplane_fix=None,
         max_nside=None,
         seeds=None,
         synalm_lmax=None,
@@ -42,7 +44,15 @@ def __init__(
             or a string (e.g. "1500 MHz") compatible with GHz.
         amplitude_modulation_temp_alm, amplitude_modulation_pol_alm: `pathlib.Path`
             Paths to the Alm expansion of the modulation maps used to rescale the small scales
-            to make them more un-uniform, they are derived from highly smoothed input emission.
+            to make them more un-uniform.
+        amplitude_modulation_beta_alm: `pathlib.Path`
+            Potentially, a different modulation map to be used for beta
+        galplane_fix: `pathlib.Path`
+            Set to None to skip the galactic plane fix in order to save some memory and
+            computing time. Used to replace the galactic emission
+            inside the GAL 070 Planck mask to a precomputed map. This is used to avoid
+            excess power in the full sky spectra due to the generated small scales being
+            too strong on the galactic plane.
         small_scale_cl, small_scale_cl_pl_index: `pathlib.Path`
             Paths to the power spectra of the small scale fluctuations for logpoltens iqu and
             the spectral index
@@ -74,7 +84,17 @@ def __init__(
                     amplitude_modulation_pol_alm,
                 ]
             ]
-            self.small_scale_cl = self.read_cl(small_scale_cl).to(u.uK_RJ ** 2)
+            if amplitude_modulation_beta_alm is not None:
+                self.modulate_alm.append(
+                    self.read_alm(amplitude_modulation_beta_alm, has_polarization=False)
+                )
+            self.small_scale_cl = self.read_cl(small_scale_cl).to(u.uK_RJ**2)
+
+        if galplane_fix is not None:
+            self.galplane_fix_map = self.read_map(galplane_fix, field=(0, 1, 2, 3))
+        else:
+            self.galplane_fix_map = None
+
         self.largescale_alm_pl_index = self.read_alm(
             largescale_alm_pl_index,
             has_polarization=False,
@@ -109,8 +129,7 @@ def draw_realization(self, synalm_lmax=None, seeds=None):
         )
 
         alm_small_scale = [
-            hp.almxfl(each, np.ones(output_lmax+1))
-            for each in alm_small_scale
+            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)
 
@@ -120,6 +139,9 @@ def draw_realization(self, synalm_lmax=None, seeds=None):
         map_small_scale[0] *= modulate_map_I
         map_small_scale[1:] *= hp.alm2map(self.modulate_alm[1].value, self.nside)
 
+        if len(self.modulate_alm) == 3:
+            modulate_map_I = hp.alm2map(self.modulate_alm[2].value, self.nside)
+
         map_small_scale += hp.alm2map(
             self.template_largescale_alm.value,
             nside=self.nside,
@@ -130,6 +152,17 @@ def draw_realization(self, synalm_lmax=None, seeds=None):
             * self.template_largescale_alm.unit
         )
 
+        if self.galplane_fix_map is not None:
+            output_IQU *= hp.ud_grade(self.galplane_fix_map[3].value, self.nside)
+            output_IQU += (
+                hp.ud_grade(
+                    self.galplane_fix_map[:3].value
+                    * (1 - self.galplane_fix_map[3].value),
+                    self.nside,
+                )
+                * self.galplane_fix_map.unit
+            )
+
         np.random.seed(seeds[1])
         output_unit = np.sqrt(1 * self.small_scale_cl_pl_index.unit).unit
         alm_small_scale = hp.synalm(
@@ -138,9 +171,7 @@ def draw_realization(self, synalm_lmax=None, seeds=None):
             new=True,
         )
 
-        alm_small_scale = hp.almxfl(
-            alm_small_scale, np.ones(output_lmax + 1)
-        )
+        alm_small_scale = hp.almxfl(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
         pl_index += (