check bugs

src/CMM/costfunc/chi2.py

@@ -76,11 +76,11 @@ def __init__(self, preset, dsim, parametric=True, full_beta_map=None):
     def _fill_A(self, x):
 
         fsub = int(
-            self.nsub * 2 / self.preset.fg.params_foregrounds["bin_mixing_matrix"]
+            self.nsub * 2 / self.preset.comp.params_foregrounds["bin_mixing_matrix"]
         )
         A = np.ones((self.nsub * 2, self.nc))
         k = 0
-        for i in range(self.preset.fg.params_foregrounds["bin_mixing_matrix"]):
+        for i in range(self.preset.comp.params_foregrounds["bin_mixing_matrix"]):
             for j in range(1, self.nc):
                 A[i * fsub : (i + 1) * fsub, j] = np.array([x[k]] * fsub)
                 k += 1
@@ -89,7 +89,7 @@ def _fill_A(self, x):
     def _get_mixingmatrix(self, nus, x):
 
         ### Compute mixing matrix
-        mixingmatrix = mm.MixingMatrix(*self.preset.fg.components_model_out)
+        mixingmatrix = mm.MixingMatrix(*self.preset.comp.components_model_out)
         return mixingmatrix.eval(nus, *x)
 
     def __call__(self, x):
@@ -121,7 +121,7 @@ def __call__(self, x):
                 )
 
                 ### Compute Planck part of the chi^2
-                mycomp = self.preset.fg.components_iter.copy()
+                mycomp = self.preset.comp.components_iter.copy()
                 # seenpix_comp = np.tile(self.preset.sky.seenpix_qubic, (mycomp.shape[0], 3, 1)).reshape(mycomp.shape)
                 mycomp[:, ~self.preset.sky.seenpix_qubic, :] = 0
                 ysim_pl = H_planck(mycomp)
@@ -189,9 +189,9 @@ def __call__(self, x):
                 )
 
                 ### Compute Planck part of the chi^2
-                # mycomp = self.preset.fg.components_iter.copy()
+                # mycomp = self.preset.comp.components_iter.copy()
                 # seenpix_comp = np.tile(self.preset.sky.seenpix_qubic, (mycomp.shape[0], 3, 1)).reshape(mycomp.shape)
-                ysim_pl = H_planck(self.preset.fg.components_iter.copy())
+                ysim_pl = H_planck(self.preset.comp.components_iter.copy())
 
                 ### Compute residuals in time domain
                 _residuals = np.r_[ysim, ysim_pl] - self.dobs
@@ -252,10 +252,10 @@ def __init__(self, preset, dsim, parametric=True, full_beta_map=None):
 
     def _fill_A(self, x):
 
-        fsub = int(self.nf / self.preset.fg.params_foregrounds["bin_mixing_matrix"])
+        fsub = int(self.nf / self.preset.comp.params_foregrounds["bin_mixing_matrix"])
         A = np.ones((self.nf, self.nc))
         k = 0
-        for i in range(self.preset.fg.params_foregrounds["bin_mixing_matrix"]):
+        for i in range(self.preset.comp.params_foregrounds["bin_mixing_matrix"]):
             for j in range(1, self.nc):
                 A[i * fsub : (i + 1) * fsub, j] = np.array([x[k]] * fsub)
                 k += 1
@@ -264,7 +264,7 @@ def _fill_A(self, x):
     def _get_mixingmatrix(self, nus, x):
 
         ### Compute mixing matrix
-        mixingmatrix = mm.MixingMatrix(*self.preset.fg.components_model_out)
+        mixingmatrix = mm.MixingMatrix(*self.preset.comp.components_model_out)
         return mixingmatrix.eval(nus, *x)
 
     def __call__(self, x):
@@ -296,7 +296,7 @@ def __call__(self, x):
                 )
 
                 ### Compute Planck part of the chi^2
-                mycomp = self.preset.fg.components_iter.copy()
+                mycomp = self.preset.comp.components_iter.copy()
                 seenpix_comp = np.tile(
                     self.preset.sky.seenpix_qubic, (mycomp.shape[0], 3, 1)
                 ).reshape(mycomp.shape)
@@ -352,9 +352,9 @@ def __call__(self, x):
                 )
 
                 ### Compute Planck part of the chi^2
-                # mycomp = self.preset.fg.components_iter.copy()
+                # mycomp = self.preset.comp.components_iter.copy()
                 # seenpix_comp = np.tile(self.preset.sky.seenpix_qubic, (mycomp.shape[0], 3, 1)).reshape(mycomp.shape)
-                ysim_pl = H_planck(self.preset.fg.components_iter.copy())
+                ysim_pl = H_planck(self.preset.comp.components_iter.copy())
 
                 ### Compute residuals in time domain
                 _residuals = np.r_[ysim, ysim_pl] - self.dobs
@@ -406,15 +406,15 @@ def __init__(self, preset, d, betamap, seenpix_wrap=None):
 
             index_num = hp.ud_grade(
                 self.preset.sky.seenpix_qubic,
-                self.preset.fg.params_foregrounds["Dust"]["nside_beta_out"],
+                self.preset.comp.params_foregrounds["Dust"]["nside_beta_out"],
             )  #
             index = np.where(index_num == True)[0]
             self._index = index
             self.seenpix_wrap = seenpix_wrap
             self.constant = False
 
     def _get_mixingmatrix(self, nus, x):
-        mixingmatrix = mm.MixingMatrix(*self.preset.fg.components_model_out)
+        mixingmatrix = mm.MixingMatrix(*self.preset.comp.components_model_out)
 
         # if self.constant:
         return mixingmatrix.eval(nus, *x)
@@ -469,9 +469,9 @@ def __call__(self, x):
             self.betamap,
             gain=self.preset.gain.gain_iter,
             fwhm=self.preset.acquisition.fwhm_mapmaking,
-            nu_co=self.preset.fg.nu_co,
+            nu_co=self.preset.comp.nu_co,
         ).operands[1]
-        tod_pl_s = H_planck(self.preset.fg.components_iter)
+        tod_pl_s = H_planck(self.preset.comp.components_iter)
 
         _r_pl = self.preset.acquisition.TOD_external - tod_pl_s
 
@@ -495,9 +495,9 @@ def __init__(self, preset, d, A_blind, icomp, seenpix_wrap=None):
         self.icomp = icomp
         self.nsub = self.preset.qubic.joint_out.qubic.Nsub
         self.fsub = int(
-            self.nsub * 2 / self.preset.fg.params_foregrounds["bin_mixing_matrix"]
+            self.nsub * 2 / self.preset.comp.params_foregrounds["bin_mixing_matrix"]
         )
-        self.nc = len(self.preset.fg.components_out)
+        self.nc = len(self.preset.comp.components_out)
 
         self.constant = True
 
@@ -522,14 +522,14 @@ def __init__(self, preset, d, A_blind, icomp, seenpix_wrap=None):
         #     self.dfg220 = self.d220[:, 1, :].copy()
         #     self.npixnf, self.nc, self.nsnd = self.d150.shape
 
-        #     index_num = hp.ud_grade(self.preset.sky.seenpix_qubic, self.preset.fg.params_foregrounds['Dust']['nside_fit'])    #
+        #     index_num = hp.ud_grade(self.preset.sky.seenpix_qubic, self.preset.comp.params_foregrounds['Dust']['nside_fit'])    #
         #     index = np.where(index_num == True)[0]
         #     self._index = index
         #     self.seenpix_wrap = seenpix_wrap
         #     self.constant = False
 
     def _get_mixingmatrix(self, x):
-        mixingmatrix = mm.MixingMatrix(self.preset.fg.components_out[self.icomp])
+        mixingmatrix = mm.MixingMatrix(self.preset.comp.components_out[self.icomp])
         if self.constant:
             return mixingmatrix.eval(self.preset.qubic.joint_out.qubic.allnus, *x)
         else:
@@ -539,7 +539,7 @@ def get_mixingmatrix_comp(self, x):
         A_comp = self._get_mixingmatrix(x)
         A_blind = self.A_blind
         print("test", A_comp.shape, A_blind.shape)
-        for ii in range(self.preset.fg.params_foregrounds["bin_mixing_matrix"]):
+        for ii in range(self.preset.comp.params_foregrounds["bin_mixing_matrix"]):
             A_blind[ii * self.fsub : (ii + 1) * self.fsub, self.icomp] = A_comp[
                 ii * self.fsub : (ii + 1) * self.fsub
             ]
@@ -601,8 +601,8 @@ def __call__(self, x):
         # H_planck = self.preset.qubic.joint_out.get_operator(self.betamap,
         #                                             gain=self.preset.gain.gain_iter,
         #                                             fwhm=self.preset.acquisition.fwhm_mapmaking,
-        #                                             nu_co=self.preset.fg.nu_co).operands[1]
-        # tod_pl_s = H_planck(self.preset.fg.components_iter)
+        #                                             nu_co=self.preset.comp.nu_co).operands[1]
+        # tod_pl_s = H_planck(self.preset.comp.components_iter)
 
         # _r_pl = self.preset.acquisition.TOD_external - tod_pl_s
         # LLH = _dot(_r.T, self.preset.acquisition.invN.operands[0](_r), self.preset.comm) + _r_pl.T @ self.preset.acquisition.invN.operands[1](_r_pl)
@@ -615,7 +615,7 @@ class Chi2Blind:
     def __init__(self, preset):
 
         self.preset = preset
-        self.nc = len(self.preset.fg.components_out)
+        self.nc = len(self.preset.comp.components_out)
         self.nf = self.preset.qubic.joint_out.qubic.nsub
         self.nsnd = (
             self.preset.qubic.joint_out.qubic.ndets
@@ -633,11 +633,11 @@ def _reshape_A(self, x):
     def _fill_A(self, x):
 
         fsub = int(
-            self.nsub * 2 / self.preset.fg.params_foregrounds["bin_mixing_matrix"]
+            self.nsub * 2 / self.preset.comp.params_foregrounds["bin_mixing_matrix"]
         )
         A = np.ones((self.nsub * 2, self.nc - 1))
         k = 0
-        for i in range(self.preset.fg.params_foregrounds["bin_mixing_matrix"]):
+        for i in range(self.preset.comp.params_foregrounds["bin_mixing_matrix"]):
             for j in range(self.nc - 1):
                 A[i * fsub : (i + 1) * fsub, j] = np.array([x[k]] * fsub)
                 k += 1
@@ -646,11 +646,11 @@ def _fill_A(self, x):
     def _reshape_A_transpose(self, x):
 
         fsub = int(
-            self.nsub * 2 / self.preset.fg.params_foregrounds["bin_mixing_matrix"]
+            self.nsub * 2 / self.preset.comp.params_foregrounds["bin_mixing_matrix"]
         )
         x_reshape = np.ones(self.nsub * 2)
 
-        for i in range(self.preset.fg.params_foregrounds["bin_mixing_matrix"]):
+        for i in range(self.preset.comp.params_foregrounds["bin_mixing_matrix"]):
             x_reshape[i * fsub : (i + 1) * fsub] = np.array([x[i]] * fsub)
         return x_reshape
 

src/CMM/params.yml

@@ -99,10 +99,10 @@ PCG:
 
     ### This section control PCG parameters
 
-    n_init_iter_pcg: 5 
-    n_iter_pcg: 5                                     # Number of PCG iterations
+    n_init_iter_pcg: 10 
+    n_iter_pcg: 10                                     # Number of PCG iterations
     tol_pcg: 1.0e-20                                # Tolerance for PCG
-    n_iter_loop: 10                                  # Number of loop (1 loop is PCG + beta fitting + gain reconstruction)
+    n_iter_loop: 100                                  # Number of loop (1 loop is PCG + beta fitting + gain reconstruction)
     ites_to_converge: 1                            # should be less than k
     tol_rms: 1.0e-10
     #fix_pixels_outside_patch: True         # Fixing pixels outside QUBIC patch

src/CMM/pipeline.py

@@ -31,7 +31,7 @@
 )
 
 # from simtools.analysis import *
-from .preset import *
+from .preset.preset import *
 
 
 class Pipeline:
@@ -128,7 +128,7 @@ def fisher_compute_sigma_r(self):
 
         # Apply Gaussian beam convolution
         C = HealpixConvolutionGaussianOperator(
-            fwhm=self.preset.acquisition.fwhm_reconstructed
+            fwhm=self.preset.acquisition.fwhm_rec
         )
         map_to_namaster = C(
             self.preset.comp.components_iter[0] - self.preset.comp.components_out[0]
@@ -140,7 +140,7 @@ def fisher_compute_sigma_r(self):
         # Compute power spectra using NaMaster
         leff, cls, _ = self.preset.sky.namaster.get_spectra(
             map_to_namaster.T,
-            beam_correction=np.rad2deg(self.preset.acquisition.fwhm_reconstructed),
+            beam_correction=np.rad2deg(self.preset.acquisition.fwhm_rec),
             pixwin_correction=False,
             verbose=False,
         )
@@ -175,14 +175,14 @@ def call_pcg(self, max_iterations, seenpix):
         initial_maps = self.preset.comp.components_iter[:, seenpix, :].copy()
 
         ### Update the precondtionner M
-        self.preset.acquisition.M = self.preset.acquisition._get_preconditioner(
+        self.preset.acquisition.M = self.preset.acquisition.get_preconditioner(
             A_qubic=self.preset.acquisition.Amm_iter[
                 : self.preset.qubic.params_qubic["nsub_out"]
             ],
             A_ext=self.preset.mixingmatrix.Amm_in[
                 self.preset.qubic.params_qubic["nsub_out"] :
             ],
-            precond=self.preset.qubic.params_qubic["preconditionner"],
+            precond=self.preset.qubic.params_qubic["preconditioner"],
         )
 
         ### Run PCG
@@ -194,7 +194,7 @@ def call_pcg(self, max_iterations, seenpix):
             maxiter = max_iterations
 
         if self.preset.tools.params["PCG"]["do_gif"]:
-            gif_folder = f"src/CMM/jobs/{self.preset.job_id}/iter/"
+            gif_folder = f"CMM/jobs/{self.preset.job_id}/iter/"
         else:
             gif_folder = None
 
@@ -226,7 +226,7 @@ def call_pcg(self, max_iterations, seenpix):
         if self.preset.tools.rank == 0:
             if self.preset.tools.params["PCG"]["do_gif"]:
                 do_gif(
-                    f"src/CMM/jobs/{self.preset.job_id}/iter/",
+                    f"CMM/jobs/{self.preset.job_id}/iter/",
                     "iter_",
                     output="animation.gif",
                 )
@@ -1179,15 +1179,15 @@ def save_data(self, step):
 
                         if step != 0:
                             os.remove(
-                                "src/CMM/"
+                                "CMM/"
                                 + self.preset.tools.params["foldername"]
                                 + "/maps/"
                                 + self.preset.tools.params["filename"]
                                 + f"_seed{str(self.preset.tools.params['CMB']['seed'])}_{str(self.preset.job_id)}_k{step-1}.pkl"
                             )
 
                     with open(
-                        "src/CMM/"
+                        "CMM/"
                         + self.preset.tools.params["foldername"]
                         + "/maps/"
                         + self.preset.tools.params["filename"]
@@ -1214,9 +1214,9 @@ def save_data(self, step):
                                 "seenpix": self.preset.sky.seenpix,
                                 "fwhm_in": self.preset.acquisition.fwhm_tod,
                                 "fwhm_out": self.preset.acquisition.fwhm_mapmaking,
-                                "fwhm_rec": self.preset.acquisition.fwhm_reconstructed,
+                                "fwhm_rec": self.preset.acquisition.fwhm_rec,
                                 #'fwhm':self.preset.acquisition.fwhm_tod,
-                                #'acquisition.fwhm_reconstructed':self.preset.acquisition.fwhm_mapmaking
+                                #'acquisition.fwhm_rec':self.preset.acquisition.fwhm_mapmaking
                             },
                             handle,
                             protocol=pickle.HIGHEST_PROTOCOL,