Merge pull request #39 from simonsobs/202208_galactic_foregrounds

New PySM 3 Galactic foregrounds simulation plan

LOGBOOK.md

@@ -2,6 +2,11 @@
 
 Date based updates from the map based simulation team
 
+## 21 March 2023 - Released MBS simulation 12
+
+Galactic, extragalactic and CMB with realistic bandpasses.
+For noise simulations, see <https://github.com/simonsobs/mnms>
+
 ## 18 Jun 2020 - Released fixed noise simulation
 
 Released new noise simulation, partly as maps on disk (currently only 1 with splits) and partly to be simulated on the fly with the `mapsims` package. See <https://github.com/simonsobs/map_based_simulations/tree/master/202006_noise> for more details.

README.md

@@ -2,7 +2,7 @@ Map based simulations for the Simons Observatory
 ================================================
 
 Map based simulations for Simons Observatory,
-maintained by the Map-Based Simulation Pipeline Working Group (MBS), [Wiki page, restricted access](http://simonsobservatory.wikidot.com/pwg:mbs)
+maintained by the Map-Based Simulation Pipeline Working Group (MBS)
 
 Github projects page for ongoing work and collaboration: https://github.com/orgs/simonsobs/projects/8
 
@@ -12,16 +12,16 @@ This repository hosts configuration files and documentation for the simulation r
 
 Maintained repositories:
 
-* [`so_pysm_models`: PySM models for Simons Observatory](https://github.com/simonsobs/so_pysm_models)
-* [`mapsims`: Map-based simulation package based on PySM + SO noise simulations](https://github.com/simonsobs/mapsims)
-* [`pysm`: Development version of PySM 3](https://github.com/healpy/pysm)
+* [`mapsims`: Map-based simulation package based on PySM + SO noise simulations](https://github.com/galsci/mapsims)
+* [`pysm`: Development version of PySM 3](https://github.com/galsci/pysm)
 
 ## Logbook
 
 Check [the Logbook](LOGBOOK.md) for the latest information about Map Based Simulations.
 
 ## Available simulations
 
+* [`mbs-s0012-20230321`](mbs-s0012-20230321/README.md): Galactic, extragalactic, CMB with realistic bandpasses.
 * [`mbs-s0011-20200618`](202006_noise/README.md): Reran of realistic noise version 3.1.1 with MSS1 hitmaps
 * [`mbs-s0010-20200221`](202002_noise/README.md): **Broken** Realistic noise version 3.1.1 with MSS1 hitmaps
 * [`mbs-s0009-20191107`](201911_lensed_cmb/README.md): 100 CMB realizations lensed with gaussian potential
@@ -33,5 +33,3 @@ Check [the Logbook](LOGBOOK.md) for the latest information about Map Based Simul
 * [`mbs-s0003-20190517`](201904_highres_foregrounds_equatorial/README.md): High resolution foregrounds in Equatorial with spectral index not varying spatially
 * [`mbs-s0002-20190325`](201903_highres_foregrounds/README.md): High resolution foregrounds with spectral index not varying spatially
 * [`mbs-s0001-20190220`](201901_gaussian_fg_lensed_cmb_realistic_noise/README.md): Gaussian foregrouds, lensed CMB, realistic noise map based simulations
-
-List of releases also available on the [`maps2params` page, restricted](http://simonsobservatory.wikidot.com/cross-group:maps2params)

mbs-s0012-20230321/PLAN.md

@@ -0,0 +1,81 @@
+Foreground simulation with realistic bandpasses and new Dust/Synchrotron models
+===============================================================================
+
+Produce full sky simulations of Galactic foregrounds using PySM 3 with the new Dust and Synchrotron models, using realistic bandpasses instead of top-hat.
+CMB and Extragalactic will be in a future release, [waiting for the new run of Websky](https://github.com/galsci/pysm/issues/120)
+
+## Foreground models
+
+The plan is to run 3 different sky models, optimistic/baseline/pessimistic, see the discussion in [this PySM 3 issue](https://github.com/galsci/pysm/issues/103#issuecomment-1081241879), for convenience here the 3 sets of models:
+
+1. `optimistic`: d9,s4,f1,a1,co1
+2. `bestguess`: d10,s5,f1,a1,co3
+3. `pessimistic`: d12,s7,f1,a2,co3
+
+Documentation reference:
+
+* `d9` `d10` GNILC based models and `d12` MKD 3D layered dust model: https://pysm3.readthedocs.io/en/latest/models.html#dust
+* Synchrotron models `s4` and `s5`: https://pysm3.readthedocs.io/en/latest/models.html#synchrotron
+* CO: https://pysm3.readthedocs.io/en/latest/models.html#co-line-emission
+* All other models are the same of PySM 2: https://pysm3.readthedocs.io/en/latest/models.html
+
+There are 2 issues:
+
+### Rotation to Equatorial
+
+In previous runs I had created Equatorial versions of all the templates to avoid having to rotate the maps because it was not supported by libsharp, for this release we can stop using libsharp and instead use all the Galactic templates and then rotate to Equatorial at the last step when we are in spherical harmonics space for smoothing.
+
+### Low resolution models
+
+* d9, d10, s4, s5 and s7 are available up to 8192
+* d12 to 2048
+* f1, a1, a2 to 512, however we have an extrapolation to 4096 by @NicolettaK, see <https://portal.nersc.gov/project/cmb/so_pysm_models_data/>
+* co? to 512 (2048 are available, but maps are smoothed at 1 deg anyway)
+
+I suggest to not use the extrapolated map by @NicolettaK, because they use a different method from the new Dust and Synchrotron models. We plan in the future to have AME and Free-free maps at 8192 with small scales produced in a similar fashion to the `d10`/`s5` models.
+
+I would run the maps in PySM at `min(max(2048, 2*Nside), available_resolution)`, then in the smoothing process, we can do a `alm2map` at the target resolution.
+Never use `hp.ud_grade`.
+
+## Bandpasses
+
+Use bandpasses from the [Instrument model](https://github.com/simonsobs/instrument_model/tree/master/instrument_hardware/modeled_bandpasses)
+
+Using the same configuration of the SO v3 sensitivity calculator:
+
+* *LF*: On-chip bandpass filters
+* *MF*: On-chip bandpass filters + Gain effects
+* *UHF*: On-chip bandpass filters + Gain effects
+
+Will be reformatted and added to instrument model in text format like: <https://github.com/simonsobs/mapsims/tree/master/mapsims/data/simonsobs_instrument_parameters_2020.06>
+
+## Beams
+
+Will be extracted from `sotodlib` into a table in txt format and added to `mapsims`, like [`simonsobs_instrument_parameters_2020.06.tbl`](https://github.com/simonsobs/mapsims/blob/master/mapsims/data/simonsobs_instrument_parameters_2020.06/simonsobs_instrument_parameters_2020.06.tbl) in the last release.
+
+## Pixelization
+
+Ideally we want to produce in 1 run only both HEALPix and CAR, once we have alms of the output map in equatorial, we can produce both final output HEALPix and CAR.
+
+Need to implement this into `mapsims`.
+
+Variable Nside: <https://github.com/simonsobs/mapsims/blob/master/mapsims/data/so_default_resolution.csv>
+
+Using `ell_max = 2.5 Nside`, write in the FITS headers the actual `ell_max` and `Nside` used for modelling.
+
+## Output files
+
+FITS files, full sky, Equatorial, `uK_CMB`, always IQU, even if component is I only to avoid broadcasting mistakes, single precision (float32), for HEALPix RING ordering.
+[Add more metadata to FITS headers](https://github.com/simonsobs/map_based_simulations/issues/38)
+
+1 map saved in FITS format for each:
+
+* sky model (optimistic, bestguess, pessimistic)
+* sky component (dust, synchrotron, free-free, CO - only if present in band, AME, their sum)
+* frequency channel (6 LAT, 6 SAT)
+* pixelization (HEALPix and CAR)
+
+Total 3 * 6 * 12 * 2 = 432 maps
+
+Naming conventions, is there any standard SO naming convention we should follow?
+Otherwise we can use the [same used for the last noise simulation](https://github.com/simonsobs/map_based_simulations/tree/master/202006_noise#available-maps)

mbs-s0012-20230321/README.md

@@ -0,0 +1,135 @@
+# Map based simulation 12: Galactic, Extragalactic, CMB with realistic bandpasses
+
+Tag: `mbs-s0012-20230321`
+
+## Updates
+
+* 2023-03-21: first release
+
+## Summary
+
+Full sky simulations for all Simons Observatory frequency channels for LAT and SAT of Galactic foregrounds using PySM 3 with the new Dust and Synchrotron models, using realistic bandpasses instead of top-hat.
+
+## Sky model
+
+This release is based on the [3 sets of recommended sky models by the Panexperiment Galactic science group](https://galsci.github.io/blog/2022/common-fiducial-sky/), in summary:
+
+* Low complexity `d9,s4,f1,a1,co1`
+* Medium complexity `d10,s5,f1,a1,co3`
+* High complexity `d12,s7,f1,a2,co3`
+
+and based on Websky for extragalactic and CMB:
+
+* `cib1,ksz1,tsz1,rg1`, `rg` stands Radio Galaxies.
+* `c3`: CMB with same Cosmological parameters used in Websky unlensed
+* `c4`: Same as `c3` but lensed by Websky
+
+Documentation reference:
+
+* `d9` `d10` GNILC based models and `d12` MKD 3D layered dust model: https://pysm3.readthedocs.io/en/latest/models.html#dust
+* Synchrotron models `s4` and `s5`: https://pysm3.readthedocs.io/en/latest/models.html#synchrotron
+* CO: https://pysm3.readthedocs.io/en/latest/models.html#co-line-emission
+* All other Galactic models are the same of PySM 2: https://pysm3.readthedocs.io/en/latest/models.html
+* For Extragalactic and CMB see [the PySM 3 documentation about Websky](https://pysm3.readthedocs.io/en/latest/websky.html#websky)
+
+## Instrument model
+
+See the `simonsobs_instrument_parameters_2023.03/` folder for test files which encode all the instrument model parameters extracted from `sotodlib` and the [notebook used to generate them](simonsobs_instrument_parameters_2023.03/extract_so_instrument_parameters.ipynb).
+
+Realistic bandpasses from the [Instrument model](https://github.com/simonsobs/instrument_model/tree/master/instrument_hardware/modeled_bandpasses), see plots in the above-mentioned notebook.
+Using the same configuration of the SO v3 sensitivity calculator:
+
+* *LF*: On-chip bandpass filters
+* *MF*: On-chip bandpass filters + Gain effects
+* *UHF*: On-chip bandpass filters + Gain effects
+
+## Available maps
+
+Maps are available both in HEALPix and in CAR (Fejer1 variant) pixelizations, generated from the same set of Alms. The resolution of the maps varies by channel, all resolutions are available in the main instrument model table.
+
+Maps are in Equatorial Coordinates, `uK_CMB` units, FITS format.
+See [`common.toml`](common.toml) for the naming convention.
+
+Each of the 17 components is available separately, see the TOML files in this repository for the configuration used to run PySM for each component.
+
+The available combination maps are 4, (see the [`combine_maps.py` script](combine_maps.py)):
+
+* `galactic_foregrounds_mediumcomplexity`
+* `galactic_foregrounds_lowcomplexity`
+* `galactic_foregrounds_highcomplexity`
+* `extragalactic`
+
+Which are meant to be used with either `cmb` for the Lensed CMB or with `cmb_unlensed`.
+
+In case you only need 1 single set of maps with all the components, you should sum `galactic_foregrounds_mediumcomplexity`, `cmb` and `extragalactic`.
+
+**Location at NERSC**, this folder on the Simons Observatory project space only includes the 4 combination maps and the 2 CMB maps (total of .75TB) due to space constraints:
+
+    /project/projectdirs/sobs/v4_sims/mbs/mbs-s0012-20230321
+
+Individual components (2.2 TB) are available on Cori scratch (accessible from Cori and from the Cori JupyterHub node, it needs membership to the `sobs` group for read access):
+
+    /global/cscratch1/sd/zonca/mbs-s0012-20230321
+
+Please [open an issue here](https://github.com/galsci/pysm/issues/new) for any data access problems.
+
+## Metadata
+
+Most useful metadata is available in the FITS header of the HEALPix maps, for example:
+
+```
+TTYPE1  = 'TEMPERATURE'                                                         
+TFORM1  = '1024E   '                                                            
+TUNIT1  = 'uK_CMB  '                                                            
+TTYPE2  = 'Q_POLARISATION'                                                      
+TFORM2  = '1024E   '                                                            
+TUNIT2  = 'uK_CMB  '                                                            
+TTYPE3  = 'U_POLARISATION'                                                      
+TFORM3  = '1024E   '                                                            
+TUNIT3  = 'uK_CMB  '                                                            
+PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation                           
+ORDERING= 'NESTED  '           / Pixel ordering scheme, either RING or NESTED   
+COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)   
+EXTNAME = 'xtension'           / name of this binary table extension            
+NSIDE   =                 1024 / Resolution parameter of HEALPIX                
+FIRSTPIX=                    0 / First pixel # (0 based)                        
+LASTPIX =             12582911 / Last pixel # (0 based)                         
+INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT                 
+OBJECT  = 'FULLSKY '           / Sky coverage, either FULLSKY or PARTIAL        
+TELESCOP= 'SAT     '                                                            
+BAND    = 'SAT_f230'                                                            
+TAG     = 'synchrotron'                                                         
+NUM     =                    0                                                  
+ELL_MAX =                 2560                                                  
+NSPLITS =                    1                                                  
+SPLIT   =                    1                                                  
+```
+
+## Model execution
+
+Simulations were run using `mapsims 2.6.0` to coordinate the execution of `PySM 3.4.0b8`.
+Given that each channel requested a different resolution, we have followed some guidelines, agreed with the Panexperiment Galactic science group:
+
+* We have 2 resolution parameters, the output Nside is the requested resolution of the output map as defined in the instrument model. The modeling Nside instead is the resolution used to run PySM, then the output of PySM is transformed to Alm, beam-smoothed, rotated to Equatorial and anti-transformed to the output Nside. No `ud_grade` operations are ever performed.
+* Evaluation of the PySM 3 models is executed at a minimum Nside 2048 or at the higher resolution available in the model. For example PySM 2 native models are executed at Nside 512, the new PySM 3 models are executed at 2048 even if we only want a Nside 128 output.
+* Evaluation is executed at 2 times the requested output Nside, unless the requested output Nside is already the maximum available. For example if we request output at Nside 2048, `d10` is executed at 4096, if we request Nside 8192, `d10` is also executed at 8192.
+* The maximum Ell is set to 2.5 times the lowest between the modeling and the output Nside, to avoid artifacts in the Spherical Harmonics transforms. Harmonics transforms are executed with [`hp.map2alm_lsq`](https://healpy.readthedocs.io/en/latest/generated/healpy.sphtfunc.map2alm_lsq.html) with 10 maximum iterations and 1e-7 target accuracy.
+
+## Verification
+
+Interactive power spectra plots for all components except CO and radio galaxies, the source notebook is available in the `verification/` folder:
+
+* [SAT TT](https://nbviewer.org/gist/zonca/3645fe8042c7d913213f3dbd647be0d5)
+* [LAT TT](https://nbviewer.org/gist/zonca/7026e5f4fd9ef304a89f1c171e43f2ce)
+* [SAT EE](https://nbviewer.org/gist/zonca/612defba6ad8d4137781661cda110bd9)
+* [LAT EE](https://nbviewer.org/gist/zonca/242527ea6bc1bc67085a7f2ae480e2df)
+* [SAT BB](https://nbviewer.org/gist/zonca/c08d6bcdd5b920b459eac77256c1d36b)
+* [LAT BB](https://nbviewer.org/gist/zonca/e9f3eebb8304583f9874cd61da95aeed)
+
+## Known issues
+
+* Websky Radio sources are available only down to 18.7 GHz, the lowest Simons Observatory channels have bandpasses to 10 Ghz, so I created a copy of 18.7 GHz and renamed it to 1.0 GHz. This is the border of the band, should not matter much.
+
+## Feedback
+
+If anything looks even just suspicious in the simulations, please do not hesitate to [open an issue here](https://github.com/galsci/pysm/issues/new) and attach a Notebook to easily reproduce the problem.

mbs-s0012-20230321/ame.toml

@@ -0,0 +1,5 @@
+tag = "ame"
+modeling_nside = 512
+
+[pysm_components]
+pysm_components_string = "a1"

mbs-s0012-20230321/ame_high.toml

@@ -0,0 +1,5 @@
+tag = "ame_high"
+modeling_nside = 512
+
+[pysm_components]
+pysm_components_string = "a2"

mbs-s0012-20230321/cib.toml

@@ -0,0 +1,5 @@
+tag = "cib"
+modeling_nside = 4096
+
+[pysm_components]
+pysm_components_string = "cib1"

mbs-s0012-20230321/cmb.toml

@@ -0,0 +1,5 @@
+tag = "cmb"
+modeling_nside = 4096
+
+[ pysm_components ]
+pysm_components_string = "c4"

mbs-s0012-20230321/cmb_unlensed.toml

@@ -0,0 +1,5 @@
+tag = "cmb_unlensed"
+modeling_nside = 4096
+
+[ pysm_components ]
+pysm_components_string = "c3"

mbs-s0012-20230321/co.toml

@@ -0,0 +1,5 @@
+tag = "co"
+modeling_nside = 2048
+
+[pysm_components]
+pysm_components_string = "co3"

mbs-s0012-20230321/co_low.toml

@@ -0,0 +1,5 @@
+tag = "co_low"
+modeling_nside = 2048
+
+[pysm_components]
+pysm_components_string = "co1"

mbs-s0012-20230321/combine_maps.py

@@ -0,0 +1,90 @@
+import os
+import healpy as hp
+import numpy as np
+from astropy.table import QTable
+from pixell import enmap
+
+# Dipole needs to be the last component, I remove dipole
+# from previous map before adding dipole
+all_combined = {
+    "galactic_foregrounds_mediumcomplexity": [
+        "dust",
+        "synchrotron",
+        "freefree",
+        "ame",
+        "co",
+    ],
+    "galactic_foregrounds_lowcomplexity": [
+        "dust_low",
+        "synchrotron_low",
+        "freefree",
+        "ame",
+        "co_low",
+    ],
+    "galactic_foregrounds_highcomplexity": [
+        "dust_high",
+        "synchrotron_high",
+        "freefree",
+        "ame_high",
+        "co",
+    ],
+    "combined_extragalactic": ["cib", "ksz", "tsz", "radio"],
+}
+
+chs = QTable.read(
+    "simonsobs_instrument_parameters_2023.03/simonsobs_instrument_parameters_2023.03.tbl",
+    format="ascii.ipac",
+)
+
+num = 0
+
+#for pixelization in ["healpix", "car"]:
+for pixelization in ["car"]:
+    for tag, components in all_combined.items():
+        for row in chs:
+            band = row["band"]
+            telescope = row["telescope"]
+            output_folder = f"output/{tag}/"
+            output_filename = (
+                output_folder
+                + f"sobs_mbs-s0012-20230321_{telescope}_mission_{band}_{tag}_{pixelization}.fits"
+            )
+            if not os.path.exists(output_filename):
+                print(20 * "*")
+                print("Starting")
+                print(output_filename)
+                for content in components:
+                    print(content)
+                    folder = f"output/{content}/"
+                    filename = (
+                        folder
+                        + f"sobs_mbs-s0012-20230321_{telescope}_mission_{band}_{content}_{pixelization}.fits"
+                    )
+                    print("Read", filename)
+                    if pixelization == "healpix":
+                        m = hp.read_map(
+                            filename, dtype=np.float64, field=(0, 1, 2), nest=True
+                        )
+                    else:
+                        m = enmap.read_map(filename)
+                    try:
+                        combined_map += m
+                    except NameError:
+                        combined_map = m
+                os.makedirs(os.path.dirname(output_filename), exist_ok=True)
+                if pixelization == "healpix":
+                    hp.write_map(
+                        output_filename,
+                        combined_map,
+                        nest=True,
+                        coord="C",
+                        column_units="uK_CMB",
+                        overwrite=True,
+                        dtype=np.float32,
+                    )
+                else:
+                    enmap.write_map(output_filename, combined_map)
+                del combined_map
+                print(20 * "*")
+                print(output_filename)
+                print(20 * "*")