Merge pull request #18 from simonsobs/cleanup_legacy

Cleanup: legacy scripts, renaming, updated README, pep8

.flake8

@@ -0,0 +1,5 @@
+[flake8]
+exclude =
+    build/,
+    legacy/,
+    setup.py

.github/environment.yml

@@ -0,0 +1,3 @@
+name: test
+dependencies:
+  python=3.8

.github/workflows/ci.yml

@@ -0,0 +1,24 @@
+name: continuous-integration
+
+on:
+  push:
+    branches:
+      - main
+      - releases/*
+  pull_request: null
+
+env:
+  CACHE_NUMBER: 0  # increase to reset cache manually
+  CONDA_ENV: .github/environment.yml
+
+jobs:
+  lint:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout SOOPERCOOL repository
+        uses: actions/checkout@v3
+
+      - name: Lint
+        uses: py-actions/flake8@v2
+        with:
+          args: "--config .flake8"

Installation_at_NERSC.md

@@ -1,4 +1,4 @@
-#### Installation of BB MASTER at `NERSC`
+#### Installation of SOOPERCOOL at `NERSC`
 
 In principle, a standard python environment with a working version of `NaMaster` will do the job. In a general environment (like `NERSC` or a different computing cluster), a safe and fast option is `micromamba`:
 
@@ -18,11 +18,11 @@ In principle, a standard python environment with a working version of `NaMaster`
   micromamba config set channel_priority strict
   ```
 
-- Install and load `bbmaster` environment (`python=3.6` proved to work with `conda`-like environments)
+- Install and load `soopercool` environment (`python=3.6` proved to work with `conda`-like environments)
 
   ```
-  micromamba create -n bbmaster python=3.6 numpy scipy ipython matplotlib healpy astropy -y
-  micromamba activate bbmaster
+  micromamba create -n soopercool python=3.6 numpy scipy ipython matplotlib healpy astropy -y
+  micromamba activate soopercool
   ```
 
   Then, install `GSL` from source (otherwise `NaMaster` won't recognize it):
@@ -46,7 +46,7 @@ In principle, a standard python environment with a working version of `NaMaster`
   micromamba install -c conda-forge namaster
   ```
 
-- Install `bbmaster` pipeline and its dependencies:
+- Install `soopercool` pipeline and its dependencies:
 
   ```
   git clone [email protected]:simonsobs/so_noise_models.git

README.md

@@ -1,26 +1,12 @@
 # Installation
 
-- [Install mbatch](https://github.com/simonsobs/mbatch#installing)
-
-- Install BBMASTER
+- Install SOOPERCOOL
 ```
-git clone [email protected]:simonsobs/BBMASTER.git
-cd BBMASTER
+git clone [email protected]:simonsobs/SOOPERCOOL.git
+cd SOOPERCOOL
 pip install -e .
 ```
 `-e` option will install local src with [development/editable mode](https://setuptools.pypa.io/en/latest/userguide/development_mode.html)
 
-# Run test
-- Prepare input files
-```
-cd data
-python get_bandpower_edges.py
-python gen_PL_sims.py
-gzip -dk mask_binary.fits.gz
-cd ..
-```
-
-- run `test_mbatch.yml`
-```
-mbatch project_name test_mbatch.yml
-```
+- Test SOOPERCOOL
+Go to the `pipeline` directory (`cd pipeline`), and run the current pipeline bash script (`./run_all.sh`).

bbmaster/deltasim.py → legacy/deltasim.py

@@ -5,7 +5,8 @@
 
 
 class DeltaBbl(object):
-    def __init__(self, nside, dsim, filt, bins, lmin=2, lmax=None, nsim_per_ell=10, seed0=1000, n_iter=0):
+    def __init__(self, nside, dsim, filt, bins, lmin=2, lmax=None,
+                 nsim_per_ell=10, seed0=1000, n_iter=0):
         if not isinstance(dsim, dict):
             raise TypeError("For now delta simulators can only be "
                             "specified through a dictionary.")
@@ -37,13 +38,14 @@ def __init__(self, nside, dsim, filt, bins, lmin=2, lmax=None, nsim_per_ell=10,
 
     def _prepare_filtering(self):
         # Match pixel resolution
-        self.filt_d['mask'] = hp.ud_grade(self.filt_d['mask'], nside_out=self.nside)
+        self.filt_d['mask'] = hp.ud_grade(self.filt_d['mask'],
+                                          nside_out=self.nside)
 
     def _gen_gaussian_alm(self, ell):
         cl = np.zeros(3*self.nside)
         cl[ell] = 1
-        # TODO: we can save time on the SHT massively, since in this case there is no
-        # sum over ell!
+        # TODO: we can save time on the SHT massively, since in
+        # this case there is no sum over ell!
         return hp.synfast(cl, self.nside)
 
     def _gen_Z2_alm(self, ell):
@@ -52,15 +54,16 @@ def _gen_Z2_alm(self, ell):
         # Generate Z2 numbers (one per m)
         # TODO: Is it clear that it's best to excite all m's
         # rather than one (or some) at a time?
-        rans = self._oosqrt2*(2*np.random.binomial(1, 0.5,
-                                                   size=2*(ell+1))-1).reshape([2,
-                                                                               ell+1])
+        rans = self._oosqrt2 *\
+            (2*np.random.binomial(1, 0.5,
+                                  size=2*(ell+1))-1).reshape([2,
+                                                              ell+1])
         # Correct m=0 (it should be real and have twice as much variance
         rans[0, 0] *= self._sqrt2
         rans[1, 0] = 0
         # Populate alms and transform to map
-        # TODO: we can save time on the SHT massively, since in this case there is no
-        # sum over ell!
+        # TODO: we can save time on the SHT massively, since in this
+        # case there is no sum over ell!
         alms = np.zeros(self.alm_ord.getsize(3*self.nside-1),
                         dtype='complex128')
         alms[idx] = rans[0]+1j*rans[1]
@@ -101,5 +104,5 @@ def get_ells(self):
         return np.arange(self.lmin, self.lmax+1)
 
     def gen_Bbl_all(self):
-        return np.array([self.gen_Bbl_at_ell(l)
-                         for l in self.get_ells()]).T
+        return np.array([self.gen_Bbl_at_ell(ell)
+                         for ell in self.get_ells()]).T

pipeline/coadder.py

@@ -1,11 +1,8 @@
 import argparse
-import healpy as hp
 import numpy as np
 import os
-from bbmaster import BBmeta
-from bbmaster.utils import *
-import pymaster as nmt
-from itertools import product
+from soopercool import BBmeta
+from soopercool.utils import theory_cls, get_noise_cls
 import matplotlib.pyplot as plt
 
 
@@ -15,134 +12,137 @@ def coadder(args):
     meta = BBmeta(args.globals)
     nmt_binning = meta.read_nmt_binning()
     n_bins = nmt_binning.get_n_bands()
-    field_pairs = ["EE", "EB", "BB"] # TODO: ADD T*, *T and BE
-    n_field_pairs = len(field_pairs)
+    field_pairs = ["EE", "EB", "BB"]  # TODO: ADD T*, *T and BE
     binary_mask = meta.read_mask("binary")
     fsky = np.mean(binary_mask)
-    
+
     if args.data:
         cl_dir = meta.cell_data_directory
     if args.sims:
         cl_dir = meta.cell_sims_directory
-    
+
     # Set the number of sims to loop over
     Nsims = meta.num_sims if args.sims else 1
 
     # Set the type of split-pair-averaged power spectra to generate.
     types = ["cross"]
-    if args.all: types.append("all") 
-    if args.auto: types.append("auto") 
-
-    for type in types:    
-        if args.plots: cells_plots = []
-
+    if args.all:
+        types.append("all")
+    if args.auto:
+        types.append("auto")
+
+    for type in types:
+        if args.plots:
+            cells_plots = []
+
         # Load split C_ells
         for id_sim in range(Nsims):
             cells_save = {}
             debug_ct = 0
-            for map_name1, map_name2 in meta.get_ps_names_list(type=type, coadd=False):
+            for map_name1, map_name2 in meta.get_ps_names_list(type=type,
+                                                               coadd=False):
                 map_set1, id_split1 = map_name1.split("__")
                 map_set2, id_split2 = map_name2.split("__")
                 n_pairs = meta.get_n_split_pairs_from_map_sets(
                     map_set1, map_set2, type=type
                 )
                 sim_label = f"_{id_sim:04d}" if Nsims > 1 else ""
                 cells_dict = np.load(
-                    f"{cl_dir}/decoupled_pcls_nobeam_{map_name1}_{map_name2}{sim_label}.npz"
+                    f"{cl_dir}/decoupled_pcls_nobeam_{map_name1}_{map_name2}{sim_label}.npz"  # noqa
                 )
-                for id_field_pair, field_pair in enumerate(field_pairs): 
+                for id_field_pair, field_pair in enumerate(field_pairs):
                     cells_label = f"{map_set1}__{map_set2}__{field_pair}"
-                    if not cells_label in cells_save:
+                    if cells_label not in cells_save:
                         cells_save[cells_label] = np.zeros(n_bins)
                     cells_save[cells_label] += cells_dict[field_pair]/n_pairs
                     # DEBUG
-                    if cells_label=="SAT1_f145__SAT1_f145__EE":
+                    if cells_label == "SAT1_f145__SAT1_f145__EE":
                         print(f"{map_name1}_{map_name2}")
                         debug_ct += 1
-                        debug_n_pairs = n_pairs
-            #print("------------------------------\nDEBUG:", type, debug_ct, debug_n_pairs)
-
+
             for cells_label in cells_save:
                 np.savez(
-                    f"{cl_dir}/decoupled_{type}_pcls_nobeam_{cells_label}{sim_label}.npz",
-                    clb=cells_save[cells_label], 
+                    f"{cl_dir}/decoupled_{type}_pcls_nobeam_{cells_label}{sim_label}.npz",  # noqa
+                    clb=cells_save[cells_label],
                     lb=nmt_binning.get_effective_ells()
                 )
-            if args.plots: cells_plots.append(cells_save)
-
+            if args.plots:
+                cells_plots.append(cells_save)
+
         if args.plots:
-            plot_dir_rel = meta.sims_directory_rel if Nsims > 1 else meta.map_directory_rel 
+            plot_dir_rel = meta.sims_directory_rel if Nsims > 1 \
+                else meta.map_directory_rel
             plot_dir = meta.plot_dir_from_output_dir(plot_dir_rel)
             el = nmt_binning.get_effective_ells()
-            
+
             for cells_label in cells_plots[0]:
                 # Load and plot theory CMB spectra
-                fig = plt.figure(figsize=(6,4))
+                plt.figure(figsize=(6, 4))
                 plt.title(cells_label)
                 map_set1, map_set2, field_pair = cells_label.split("__")
-                el_th, cl_th_dict = theory_cls(meta.cosmology, 
+                el_th, cl_th_dict = theory_cls(meta.cosmology,
                                                meta.lmax, lmin=meta.lmin)
                 cl_th = cl_th_dict[field_pair]
-                plt.plot(el_th, 
+                plt.plot(el_th,
                          el_th*(el_th+1)/2./np.pi*(cl_th),
-                         lw=1, ls='--', c='darkred', 
-                         label=f"theory (noiseless)")
-                
-                if type=="auto" or type=="all":
+                         lw=1, ls='--', c='darkred',
+                         label="theory (noiseless)")
+
+                if type == "auto" or type == "all":
                     # Plot noisy theory with beam-deconvolved noise
                     _, nl_th_beam_dict = get_noise_cls(
-                        fsky, meta.lmax, lmin=meta.lmin, 
-                        sensitivity_mode='baseline', 
+                        fsky, meta.lmax, lmin=meta.lmin,
+                        sensitivity_mode='baseline',
                         oof_mode='optimistic',
                         is_beam_deconvolved=True
                     )
                     freq_tag = meta.freq_tag_from_map_set(map_set1)
                     n_splits = meta.n_splits_from_map_set(map_set1)
                     nl_th_beam = nl_th_beam_dict["P"][freq_tag]*float(n_splits)
-                    plt.plot(el_th, 
+                    plt.plot(el_th,
                              el_th*(el_th+1)/2./np.pi*(cl_th+nl_th_beam),
-                             lw=1, ls='-.', c='k', 
-                             label=f"theory (noisy)")
-                
+                             lw=1, ls='-.', c='k',
+                             label="theory (noisy)")
+
                 # Plot split-pair-averaged power spectra
                 cells_sum = 0
                 for id_sim in range(Nsims):
                     cells = cells_plots[id_sim][cells_label]
                     cells_sum += cells
                     if Nsims > 1:
-                        plt.plot(el, el*(el+1)/2./np.pi*cells, 
+                        plt.plot(el, el*(el+1)/2./np.pi*cells,
                                  lw=0.3, c='tab:red', alpha=0.5)
-                plt.plot(el, el*(el+1)/2./np.pi*cells_sum/Nsims, 
+                plt.plot(el, el*(el+1)/2./np.pi*cells_sum/Nsims,
                          lw=1, c='blue', label=cells_label)
 
                 plt.xscale('log')
                 plt.yscale('log')
                 plt.ylabel(r'$\ell(\ell+1)\,C_\ell/2\pi$', fontsize=14)
                 plt.xlabel(r'$\ell$', fontsize=14)
-                plt.legend(bbox_to_anchor=(1,1))
+                plt.legend(bbox_to_anchor=(1, 1))
                 plt.savefig(
                     os.path.join(
-                        plot_dir, 
+                        plot_dir,
                         f"cells_{type}_{cells_label}.png"
                     ), bbox_inches='tight'
                 )
 
+
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(description="Pseudo-Cl calculator")
     parser.add_argument("--globals", type=str, help="Path to the yaml file")
-    parser.add_argument("--auto", action="store_true", 
+    parser.add_argument("--auto", action="store_true",
                         help="Save auto-split power spectra.")
-    parser.add_argument("--all", action="store_true", 
-                        help="Save coadded auto- and cross-split power spectra.")
-    parser.add_argument("--plots", action="store_true", 
+    parser.add_argument(
+        "--all", action="store_true",
+        help="Save coadded auto- and cross-split power spectra.")
+    parser.add_argument("--plots", action="store_true",
                         help="Plot the generated power spectra if True.")
 
     mode = parser.add_mutually_exclusive_group()
     mode.add_argument("--data", action="store_true")
     mode.add_argument("--sims", action="store_true")
-
-
-
+
     args = parser.parse_args()
 
-    coadder(args)
+    coadder(args)