General refactoring and some SO specific changes

params/Planck_TODs/param_PlanckTODs_CG.yml

@@ -54,15 +54,25 @@ general:
   chain_maps_interval: 5
 
   output_paths:
-    stats: "../cmdr4_stats/PlanckTODs_CG/"
-    chains: "../cmdr4_chains/PlanckTODs_CG/"
+    stats: "../cmdr4_stats/PlanckTODs_CG_samplegroups/"
+    chains: "../cmdr4_chains/PlanckTODs_CG_samplegroups/"
 
   logging:
     console:
       level: 'info'
     file:
       level: 'info'
-      filename: '../PlanckTODs_CG.log'
+      filename: '../PlanckTODs_CG_samplegroups.log'
+
+
+
+sampling_groups_compsep:
+  sample_amps_perpix:
+    enabled: true
+    sample_class: "amplitude_sampler_CG"
+    comps: [CMB, Synchrotron, ThermalDust, FreeFree]
+    bands: [Planck30GHz, Haslam408MHz, WMAPKa, WMAPQ1, WMAPQ2, WMAPV1, WMAPV2, Planck857GHz, Planck353GHz]
+
 
 components:
   CMB:
@@ -331,11 +341,11 @@ experiments:
         tod_files_prefix: "LFI_030_"
         num_MPI_tasks: 64 #24
         eval_nside: 512
-        data_nside: 512
         fwhm: 32.4 # [arcmin]
         freq: 28.4 # [GHz]
         fsamp: 32.5079365079365 # [Hz]
         polarization: "IQU"
+        noise_psd: "NoisePSDOof"
         processing_mask: "/mn/stornext/d5/data/duncanwa/WMAP/data/mask_proc_030_res_v5.fits"
         filelist: "/mn/stornext/d16/cmbco/bp/mathew/test/filelist_30.txt"
         detectors:
@@ -360,10 +370,8 @@ experiments:
         tod_files_prefix: "LFI_044_"
         num_MPI_tasks: 40
         eval_nside: 512
-        data_nside: 512
         fwhm: 27.1 # [arcmin]
         freq: 44.1 # [GHz]
-        fsamp: 46.5454545454545
         polarization: "IQU"
         processing_mask: "/mn/stornext/d5/data/duncanwa/WMAP/data/mask_proc_044_res_v5.fits"
         filelist: "/mn/stornext/d16/cmbco/bp/mathew/test/filelist_44.txt"
@@ -397,14 +405,11 @@ experiments:
         tod_files_prefix: "LFI_070_"
         num_MPI_tasks: 64
         eval_nside: 1024
-        data_nside: 1024
         fwhm: 13.3 # [arcmin]
         freq: 70.1 # [GHz]
-        fsamp: 78.7692307692308
         polarization: "IQU"
         processing_mask: "/mn/stornext/d5/data/duncanwa/WMAP/data/mask_proc_070_res_v5.fits"
         filelist: "/mn/stornext/d16/cmbco/bp/mathew/test/filelist_70.txt"
-        detectors:
         detectors:
           18M:
             name: "18M"

params/Planck_TODs/param_PlanckTODs_perpix.yml

@@ -21,7 +21,7 @@ general:
   sample_abs_gain_from_iter_num: 2
   sample_rel_gain: True
   sample_rel_gain_from_iter_num: 2
-  sample_temporal_gain: False
+  sample_temporal_gain: True
   sample_temporal_gain_from_iter_num: 2
 
   perform_compsep: true
@@ -44,15 +44,22 @@ general:
   init_from_chain: false
 
   output_paths:
-    chains: "../cmdr4_chains/PlanckTODs_perpix_newclass_v5_abs+rel/"
-    stats: "../cmdr4_stats/PlanckTODs_perpix_newclass_v5_abs+rel/"
+    chains: "../cmdr4_chains/PlanckTODs_perpix_sampgroups/"
+    stats: "../cmdr4_stats/PlanckTODs_perpix_sampgroups/"
 
   logging:
     console:
       level: 'debug'
     file:
       level: 'info'
-      filename: '../cmdr4_logs/output_PlanckTODs_perpix_newclass_v1.log'
+      filename: '../cmdr4_logs/output_PlanckTODs_perpix_sampgroups'
+
+sampling_groups_compsep:
+  sample_amps_perpix:
+    enabled: true
+    sample_class: "amplitude_sampler_perpix"
+    comps: [CMB, Synchrotron, ThermalDust, FreeFree]
+    bands: [Planck30GHz, Haslam408MHz, WMAPKa, WMAPQ1, WMAPQ2, WMAPV1, WMAPV2, Planck857GHz, Planck353GHz]
 
 
 components:

params/param_AkariTODs_CG.yml

@@ -28,14 +28,32 @@ general:
   noise_fknee: 0.1
   noise_alpha: -0.9
 
-  sample_corr_noise: True
-  sample_corr_noise_from_iter_num: 8
-  sample_abs_gain: True
-  sample_abs_gain_from_iter_num: 2
-  sample_rel_gain: True
-  sample_rel_gain_from_iter_num: 4
-  sample_temporal_gain: True
-  sample_temporal_gain_from_iter_num: 4
+  corr_noise:
+    sample_corr_noise: True
+    sample_corr_noise_from_iter_num: 8
+    sample_noise_params: True     # Resample noise-model params (needs sample_corr_noise).
+    sample_sigma0: True            # Re-estimate sigma0 from the n_corr-subtracted residual.
+    CG_err_tol: 1.0e-4             # Correlated-noise CG convergence tolerance.
+    CG_max_iter: 60                # Max CG iterations; 0 -> skip CG, use stationary fallback.
+    nomono: False                 # Project the per-scan monopole out of n_corr (Fortran "nomono").
+    onlymono: False               # Model the correlated noise as only the per-scan offset.
+    sigma0_decimation: 1          # Block-average factor for the sigma0 estimator.
+    psd_fit_nu_min: 0.0           # Min frequency (Hz) for noise-PSD parameter fitting.
+    psd_fit_nu_max: 2.0           # Max frequency (Hz) for noise-PSD parameter fitting.
+    psd_bin: False                # Fit every Fourier mode (False) or mode-weighted binned PSD.
+  gain_calib_downsample_time: 1.0   # [s] Averaging chunk duration for gain-calibration TODs (0 = off).
+  abs_gain:                           # Absolute (band-global) gain g0.
+    sample: true
+    sample_from_iter_num: 2
+    calibrate_against: orbital_dipole # "orbital_dipole" | "full_sky" | "sky"; per-band overridable.
+  rel_gain:                           # Per-detector relative gain.
+    sample: true
+    sample_from_iter_num: 4
+    calibrate_against: full_sky       # "orbital_dipole" | "full_sky" | "sky"; per-band overridable.
+  temporal_gain:                      # Per-scan temporal gain variations.
+    sample: true
+    sample_from_iter_num: 4
+    calibrate_against: full_sky       # "orbital_dipole" | "full_sky" | "sky"; per-band overridable.
 
   tod_processing_steps:
     Akari:
@@ -69,6 +87,7 @@ general:
   write_orb_dipole_maps_to_chain: True
   write_corr_noise_maps_to_chain: True
   write_sky_model_maps_to_chain: True
+  write_ncorr_tods_to_chain: False   # DEBUG ONLY: full per-sample n_corr TODs (very large data).
 
   output_paths:
     stats: "../cmdr4_stats/PlanckTODs_CG/"
@@ -107,7 +126,8 @@ components:
       smoothing_prior_amplitude: 2.0e-3  # uK_RJ^2
       beta: 1.56
       T: 20
-      nu0: 545
+      nu_ref: [857, 353]  # [I, QU]
+      units: "uK_RJ"
       lmax: 1024
       smoothing_scale: 0.0
       spatially_varying_MM: False
@@ -122,7 +142,8 @@ components:
       smoothing_prior_FWHM: 60.0 # arcmin
       smoothing_prior_amplitude: 2.0e-2  # uK_RJ^2
       beta: -3.1
-      nu0: 30
+      nu_ref: [0.408, 30]  # [I, QU]
+      units: "uK_RJ"
       lmax: 512
       smoothing_scale: 0.0
       spatially_varying_MM: False
@@ -137,7 +158,8 @@ components:
       smoothing_prior_FWHM: 30.0 # arcmin
       smoothing_prior_amplitude: 1.0e+9  # uK_RJ^2
       T: 7000
-      nu0: 0.408 #40
+      nu_ref: 40
+      units: "uK_RJ"
       lmax: 512
       smoothing_scale: 0.0
       spatially_varying_MM: False
@@ -443,12 +465,14 @@ experiments:
     pix_is_compressed: true
     data_type: "TOD"
     mapmaker: "CG"
-    flag_bitmask: 0b1111111111111111111111111111101
+    bad_data_bitmask: 0b1111111111111111111111111111101
     #bad_PIDs_path: "/mn/stornext/d23/cmbco/jonas/c4_testing/Commander4/badPIDs.npy"
     Fourier_times_path: "/mn/stornext/d23/cmbco/jonas/c4_testing/analysis/FFT_times_0_399999.npy"
     num_scans: 11395
     bands:
       AkariN60:
+        abs_gain:
+          calibrate_against: full_sky   # Preserve full-sky abs-gain calibration (legacy nu>380 rule).
         enabled: true
         tod_files_prefix: "AKARI_N60_v1_"
         polarization: "I"

pyproject.toml

@@ -29,6 +29,8 @@ dependencies = [
     "healpy",
     "tqdm",
     "camb",
+    "numpy-quaternion",
+    "qpoint",
     ]  # Add new Python runtime dependencies here.
 
 [project.optional-dependencies]

src/commander4/cli.py

@@ -64,9 +64,11 @@ def run_commander4(params: Bunch, params_dict: dict):
     # as Numpy will not respect a change in thread count after it has been loaded.
     import numpy as np
     import commander4.output.log as log
-    from commander4.tod_processing import process_tod, init_tod_processing, get_initial_sky
-    from commander4.compsep_processing import process_compsep, init_compsep_processing
-    from commander4.communication import receive_tod, send_tod, receive_compsep, send_compsep
+    from commander4.tod_processing import process_tod, init_tod_processing
+    from commander4.compsep_processing import process_compsep, init_compsep_processing,\
+        get_initial_sky_model
+    from commander4.communication import receive_tod, send_tod, receive_compsep, send_compsep,\
+        get_local_initial_sky
 
     # Unique seed per worldrank. Hash used for slightly improved entropy. Modulus because seed needs
     # to be 32 bit. Numpy recommends instead carrying around an instance of 'np.random.default_rng',
@@ -103,21 +105,32 @@ def run_commander4(params: Bunch, params_dict: dict):
             mpi_info.world.compsep_band_masters, root=mpi_info.world.compsep_master)
         mpi_info['world']['compsep_band_masters'] = world_compsep_band_masters_dict
 
-    ###### Sending empty data back and forth ######
+    ###### Exchanging the initial sky model ######
+    # Component separation is active iff CompSep ranks were allocated (compsep_master is then a
+    # valid world rank). This single flag replaces the old `perform_compsep` parameter.
+    compsep_active = mpi_info.world.compsep_master is not None
     curr_tod_output = None
     if mpi_info.world.color == 0:
-        # Chain #1 do TOD processing, resulting in maps_chain1 (we start with a fake output of
-        # component separation, containing a completely empty sky).
-        compsep_output_black = get_initial_sky(experiment_data)
+        # The initial sky model is built from each component's init_from / init_chain_path (else
+        # zeros). If CompSep ranks exist they build and send it (as for every later iteration);
+        # otherwise we build it locally so a sensible fixed sky is available with no CompSep ranks.
+        if compsep_active:
+            curr_compsep_output = receive_compsep(mpi_info, experiment_data, my_band_tod_id,
+                                                  mpi_info.world.compsep_band_masters)
+        else:
+            curr_compsep_output = get_local_initial_sky(mpi_info, experiment_data, params)
 
         curr_tod_output, tod_samples = process_tod(mpi_info, experiment_data,
                                                         tod_samples_chain1,
-                                                        compsep_output_black, params, 1, 1)
-        if params.general.perform_compsep:
+                                                        curr_compsep_output, params, 1, 1)
+        if compsep_active:
             send_tod(mpi_info, curr_tod_output, my_band_tod_id, mpi_info.world.compsep_band_masters)
-        curr_compsep_output = compsep_output_black
 
     elif mpi_info.world.color == 1:
+        # Send the initial sky model to TOD before receiving the first TOD output, mirroring the
+        # process_compsep -> send_compsep -> receive_tod order used inside the main loop.
+        send_compsep(mpi_info, my_band_compsep_id, get_initial_sky_model(components),
+                     mpi_info.world.tod_band_masters)
         curr_tod_output = receive_tod(mpi_info, mpi_info.world.tod_band_masters, my_band,
                                       my_band_compsep_id, curr_tod_output)
 
@@ -150,7 +163,7 @@ def run_commander4(params: Bunch, params_dict: dict):
                 logger.info(f"TOD: Rank {mpi_info.tod.rank} finished chain {chain_num}, iter "\
                             f"{iter_num} in {time.time()-t0:.2f}s. Receiving compsep results.")
             t0 = time.time()
-            if params.general.perform_compsep:
+            if compsep_active:
                 curr_compsep_output = receive_compsep(mpi_info, experiment_data,
                                                     my_band_tod_id,
                                                     mpi_info.world.compsep_band_masters)
@@ -159,7 +172,7 @@ def run_commander4(params: Bunch, params_dict: dict):
                             f"results for chain {chain_num}, iter {iter_num} "\
                             f"(time spent waiting+receiving = "\
                             f"{time.time()-t0:.1f}s).")
-            if params.general.perform_compsep:
+            if compsep_active:
                 send_tod(mpi_info, curr_tod_output, my_band_tod_id,
                         mpi_info.world.compsep_band_masters)
             if mpi_info.tod.is_master:
@@ -191,7 +204,7 @@ def run_commander4(params: Bunch, params_dict: dict):
                         f" chain {chain_num}, iter {iter_num} (time spent waiting+receiving = "\
                         f"{time.time()-t0:.1f}s).")
     # stop compsep machinery
-    if mpi_info.world.is_master:
+    if mpi_info.world.is_master and compsep_active:
         logger.info("TOD: sending STOP signal to compsep")
         mpi_info.world.comm.send(True, dest=mpi_info.world.compsep_master)