From 07a22219f9503d02ac306f6087ac030b2e10d763 Mon Sep 17 00:00:00 2001 From: TomLaclavere Date: Wed, 11 Sep 2024 17:20:54 +0200 Subject: [PATCH 1/2] WIP - correct issues with nsub in & out --- src/CMM/preset/preset_acquisition.py | 2 +- src/FMM/main.py | 6 -- src/FMM/params.yml | 20 +++--- src/FMM/pipeline.py | 92 +++++++++++++++------------- src/lib/Qcg.py | 2 + src/run.py | 2 +- 6 files changed, 65 insertions(+), 59 deletions(-) diff --git a/src/CMM/preset/preset_acquisition.py b/src/CMM/preset/preset_acquisition.py index 37ac28f..d79de64 100644 --- a/src/CMM/preset/preset_acquisition.py +++ b/src/CMM/preset/preset_acquisition.py @@ -3,7 +3,7 @@ from lib.Qacquisition import * from lib.Qnoise import * -import fgb.component_model as c +import lib.Qcomponent_model as c class PresetAcquisition: diff --git a/src/FMM/main.py b/src/FMM/main.py index 30fee95..54d4ea6 100644 --- a/src/FMM/main.py +++ b/src/FMM/main.py @@ -1,9 +1,3 @@ -import Qmixing_matrix - -print(Qmixing_matrix) -sto - - import sys from pipeline import * diff --git a/src/FMM/params.yml b/src/FMM/params.yml index 0f19c72..6fa178e 100644 --- a/src/FMM/params.yml +++ b/src/FMM/params.yml @@ -1,4 +1,4 @@ -path_out: 'test_to_remove/' +path_out: 'test_FastSim/' datafilename: 'MC' CMB: @@ -17,11 +17,11 @@ QUBIC: ### This section control the parameters related to QUBIC instrument: 'DB' # Instrumental design (DB or UWB) - npointings: 6000 # Number of time samples - nsub_in: 12 # Number of sub-acquisitions - nsub_out: 12 - nrec: 6 # Number of reconstructed frequency maps - convolution_in: False # Angular resolution in the input TOD + npointings: 100 # Number of time samples + nsub_in: 10 # Number of sub-acquisitions + nsub_out: 10 + nrec: 2 # Number of reconstructed frequency maps + convolution_in: True # Angular resolution in the input TOD convolution_out: False # Angular resolution in the reconstruction bandpass_correction: True # Correction to keep the edges of the integration band constant NOISE: @@ -46,11 +46,11 @@ PLANCK: ### This section define which external data are we using - external_data: True + external_data: False weight_planck: 0 # Weight of Planck data within the QUBIC patch level_noise_planck: 1 # Noise level for Planck data bandwidth_planck: 0.2 # Multiplicative factor to define the bandwidth of Planck's data - nsub_planck: 100 # Number of sub acquisitions for Planck's data + nsub_planck: 10 # Number of sub acquisitions for Planck's data Pipeline: mapmaking: True # Run MapMaking Pipeline @@ -60,9 +60,9 @@ PCG: ### This section control PCG parameters - n_iter_pcg: 50 # Number of PCG iterations + n_iter_pcg: 100 # Number of PCG iterations tol_pcg: 1.0e-12 # Tolerance for PCG - preconditioner: True + preconditioner: False gif: True resolution_plot: 12 fwhm_plot: 0.005 diff --git a/src/FMM/pipeline.py b/src/FMM/pipeline.py index 3b2a2d1..e977b6f 100644 --- a/src/FMM/pipeline.py +++ b/src/FMM/pipeline.py @@ -67,7 +67,8 @@ def __init__(self, comm, file, params): self.params = params ### Fsub - self.fsub = int(self.params["QUBIC"]["nsub_out"] / self.params["QUBIC"]["nrec"]) + self.fsub_in = int(self.params['QUBIC']['nsub_in'] / self.params['QUBIC']['nrec']) + self.fsub_out = int(self.params['QUBIC']['nsub_out'] / self.params['QUBIC']['nrec']) self.file = file self.plot_folder = "FMM/" + self.params["path_out"] + "png/" @@ -103,32 +104,25 @@ def __init__(self, comm, file, params): self.dict_out = self.get_dict(key="out") self.skyconfig = self.get_sky_config() - ### Joint acquisition for TOD making - self.joint_tod = JointAcquisitionFrequencyMapMaking( - self.dict_in, - self.params["QUBIC"]["instrument"], - self.params["QUBIC"]["nsub_in"], - self.params["QUBIC"]["nsub_in"], - H=None, - ) - ### Joint acquisition - if self.params["QUBIC"]["nsub_in"] == self.params["QUBIC"]["nsub_out"]: - H = self.joint_tod.qubic.H - else: - H = None - - self.joint = JointAcquisitionFrequencyMapMaking( - self.dict_out, - self.params["QUBIC"]["instrument"], - self.params["QUBIC"]["nrec"], - self.params["QUBIC"]["nsub_out"], - H=H, - ) - + self.joint = JointAcquisitionFrequencyMapMaking(self.dict_out, self.params['QUBIC']['instrument'], + self.params['QUBIC']['nrec'], + self.params['QUBIC']['nsub_out']) self.planck_acquisition143 = PlanckAcquisition(143, self.joint.qubic.scene) self.planck_acquisition217 = PlanckAcquisition(217, self.joint.qubic.scene) - self.nus_Q = self.get_averaged_nus() + self.nus_Q = self._get_averaged_nus() + + ### Joint acquisition for TOD making + if self.params['QUBIC']['nsub_in'] == self.params['QUBIC']['nsub_out']: + self.joint_tod = JointAcquisitionFrequencyMapMaking(self.dict_in, self.params['QUBIC']['instrument'], + self.params['QUBIC']['nsub_in'], + self.params['QUBIC']['nsub_in'], + H=self.joint.qubic.H) + else: + self.joint_tod = JointAcquisitionFrequencyMapMaking(self.dict_in, self.params['QUBIC']['instrument'], + self.params['QUBIC']['nsub_in'], + self.params['QUBIC']['nsub_in'], + H=None) ### Coverage map self.coverage = self.joint.qubic.subacqs[0].get_coverage() @@ -148,7 +142,7 @@ def __init__(self, comm, file, params): self.external_timeline = ExternalData2Timeline( self.skyconfig, - self.joint.qubic.allnus, + self.joint_tod.qubic.allnus, self.params["QUBIC"]["nrec"], nside=self.params["SKY"]["nside"], corrected_bandpass=self.params["QUBIC"]["bandpass_correction"], @@ -223,6 +217,20 @@ def get_components_fgb(self): dict_comps += [Synchrotron(nu0=150, beta_pl=-3)] return dict_comps + + def _get_averaged_nus(self): + + """ + + Method to average QUBIC frequencies. + + """ + + nus_eff = [] + for i in range(self.params['QUBIC']['nrec']): + nus_eff += [np.mean(self.joint.qubic.allnus[i*self.fsub_out:(i+1)*self.fsub_out])] + + return np.array(nus_eff) def get_random_value(self): """Random value @@ -278,13 +286,11 @@ def get_averaged_nus(self): """ - nus_ave = [] - for i in range(self.params["QUBIC"]["nrec"]): - nus_ave += [ - np.mean(self.joint.qubic.allnus[i * self.fsub : (i + 1) * self.fsub]) - ] + nus_eff = [] + for i in range(self.params['QUBIC']['nrec']): + nus_eff += [np.mean(self.joint.qubic.allnus[i*self.fsub_out:(i+1)*self.fsub_out])] - return np.array(nus_ave) + return np.array(nus_eff) def get_sky_config(self): """Sky configuration. @@ -491,12 +497,12 @@ def get_convolution(self): for irec in range(self.params["QUBIC"]["nrec"]): numerator_fwhm, denominator_fwhm = 0, 0 numerator_nus, denominator_nus = 0, 0 - for jsub in range(irec * self.fsub_out, (irec + 1) * self.fsub_out): + for jsub in range(irec * self.fsub_in, (irec + 1) * self.fsub_in): # Compute the expected reconstructed resolution for sub-acquisition numerator_fwhm += ( scalar_acquisition_operators[jsub] * weight_factor[jsub] - * self.fwhm_in[jsub] + * fwhm_in[jsub] ) denominator_fwhm += ( scalar_acquisition_operators[jsub] * weight_factor[jsub] @@ -544,7 +550,7 @@ def get_input_map(self): for i in range(self.params["QUBIC"]["nrec"]): m_nu_in[i] = np.mean( - self.external_timeline.m_nu[i * self.fsub : (i + 1) * self.fsub], axis=0 + self.external_timeline.m_nu[i * self.fsub_out : (i + 1) * self.fsub_out], axis=0 ) return m_nu_in @@ -583,7 +589,7 @@ def get_tod(self, noise=False): if self.params["QUBIC"]["convolution_in"]: C = HealpixConvolutionGaussianOperator( fwhm=np.min( - self.fwhm_in[irec * self.fsub : (irec + 1) * self.fsub] + self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] ) ) @@ -600,7 +606,7 @@ def get_tod(self, noise=False): if self.params["QUBIC"]["convolution_in"]: C = HealpixConvolutionGaussianOperator( fwhm=np.min( - self.fwhm_in[irec * self.fsub : (irec + 1) * self.fsub] + self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] ) ) @@ -643,7 +649,7 @@ def get_tod(self, noise=False): else: - sh_q = self.joint.qubic.ndets * self.joint.qubic.nsamples + sh_q = self.joint_tod.qubic.ndets * self.joint_tod.qubic.nsamples TOD_QUBIC = ( self.H_in_qubic(factor * self.external_timeline.m_nu).ravel() + self.noiseq @@ -666,7 +672,7 @@ def get_tod(self, noise=False): if self.params["QUBIC"]["convolution_in"]: C = HealpixConvolutionGaussianOperator( fwhm=np.min( - self.fwhm_in[irec * self.fsub : (irec + 1) * self.fsub] + self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] ) ) @@ -687,7 +693,7 @@ def get_tod(self, noise=False): if self.params["QUBIC"]["convolution_in"]: C = HealpixConvolutionGaussianOperator( fwhm=np.min( - self.fwhm_in[irec * self.fsub : (irec + 1) * self.fsub] + self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] ) ) @@ -759,8 +765,8 @@ def get_preconditioner(self): ) for irec in range(self.params["QUBIC"]["nrec"]): - imin = irec * self.fsub - imax = (irec + 1) * self.fsub + imin = irec * self.fsub_out + imax = (irec + 1) * self.fsub_out for istk in range(3): conditioner[irec, self.seenpix, istk] = 1 / ( np.sum(approx_hth[imin:imax, self.seenpix, 0], axis=0) @@ -800,6 +806,10 @@ def pcg(self, d, x0, seenpix): A = self.H_out.T * self.invN * self.H_out x_planck = self.m_nu_in * (1 - seenpix[None, :, None]) + print(self.H_out_all_pix.shapein) + print(x_planck.shape) + print(d.shape) + print(self.invN) b = self.H_out.T * self.invN * (d - self.H_out_all_pix(x_planck)) ### Preconditionning M = self.get_preconditioner() diff --git a/src/lib/Qcg.py b/src/lib/Qcg.py index 75812f0..42d8144 100644 --- a/src/lib/Qcg.py +++ b/src/lib/Qcg.py @@ -180,6 +180,8 @@ def iteration(self): self.A(self.d, self.q) alpha = self.delta / self.dot(self.d, self.q) self.x += alpha * self.d + print(self.x.shape) + print(self.input[:, self.seenpix, :].shape) _r = self.x - self.input[:, self.seenpix, :] self.rms = np.std(_r, axis=1) if self.gif is not None: diff --git a/src/run.py b/src/run.py index 6cab77e..b2e0cd1 100644 --- a/src/run.py +++ b/src/run.py @@ -11,7 +11,7 @@ if __name__ == "__main__": if simu == 'FMM': - + try: file = str(sys.argv[1]) except IndexError: From 0804c15b924bdcac7342be682989787c13895bec Mon Sep 17 00:00:00 2001 From: TomLaclavere Date: Fri, 13 Sep 2024 11:03:24 +0200 Subject: [PATCH 2/2] Correct issue for FMM --- src/FMM/params.yml | 10 +- src/FMM/pipeline.py | 289 +++++++++++++++++++++------------------- src/lib/Qcg.py | 26 ++-- src/lib/Qmap_plotter.py | 2 +- 4 files changed, 171 insertions(+), 156 deletions(-) diff --git a/src/FMM/params.yml b/src/FMM/params.yml index 6fa178e..6a00fa2 100644 --- a/src/FMM/params.yml +++ b/src/FMM/params.yml @@ -17,9 +17,9 @@ QUBIC: ### This section control the parameters related to QUBIC instrument: 'DB' # Instrumental design (DB or UWB) - npointings: 100 # Number of time samples - nsub_in: 10 # Number of sub-acquisitions - nsub_out: 10 + npointings: 1000 # Number of time samples + nsub_in: 8 # Number of sub-acquisitions + nsub_out: 4 nrec: 2 # Number of reconstructed frequency maps convolution_in: True # Angular resolution in the input TOD convolution_out: False # Angular resolution in the reconstruction @@ -37,7 +37,7 @@ SKY: ### This section control the reconstructed sky - nside: 256 # Nside of components + nside: 32 # Nside of components coverage_cut: 0.2 # Define a cut where the coverage is too low RA_center: 0 DEC_center: -57 @@ -62,7 +62,7 @@ PCG: n_iter_pcg: 100 # Number of PCG iterations tol_pcg: 1.0e-12 # Tolerance for PCG - preconditioner: False + preconditioner: True gif: True resolution_plot: 12 fwhm_plot: 0.005 diff --git a/src/FMM/pipeline.py b/src/FMM/pipeline.py index e977b6f..162c683 100644 --- a/src/FMM/pipeline.py +++ b/src/FMM/pipeline.py @@ -23,6 +23,7 @@ from .model.planck_timeline import * +''' def save_pkl(name, d): """ Function that create a file called "name" to save the dictionary "d" using Pickle package @@ -36,7 +37,7 @@ def save_pkl(name, d): """ with open(name, "wb") as handle: pickle.dump(d, handle, protocol=pickle.HIGHEST_PROTOCOL) - +''' __all__ = ["PipelineFrequencyMapMaking", "PipelineEnd2End"] @@ -67,8 +68,8 @@ def __init__(self, comm, file, params): self.params = params ### Fsub - self.fsub_in = int(self.params['QUBIC']['nsub_in'] / self.params['QUBIC']['nrec']) - self.fsub_out = int(self.params['QUBIC']['nsub_out'] / self.params['QUBIC']['nrec']) + self.fsub_in = int(self.params["QUBIC"]["nsub_in"] / self.params["QUBIC"]["nrec"]) + self.fsub_out = int(self.params["QUBIC"]["nsub_out"] / self.params["QUBIC"]["nrec"]) self.file = file self.plot_folder = "FMM/" + self.params["path_out"] + "png/" @@ -104,25 +105,32 @@ def __init__(self, comm, file, params): self.dict_out = self.get_dict(key="out") self.skyconfig = self.get_sky_config() + ### Joint acquisition for TOD making + self.joint_tod = JointAcquisitionFrequencyMapMaking( + self.dict_in, + self.params["QUBIC"]["instrument"], + self.params["QUBIC"]["nsub_in"], + self.params["QUBIC"]["nsub_in"], + H=None, + ) + ### Joint acquisition - self.joint = JointAcquisitionFrequencyMapMaking(self.dict_out, self.params['QUBIC']['instrument'], - self.params['QUBIC']['nrec'], - self.params['QUBIC']['nsub_out']) + if self.params["QUBIC"]["nsub_in"] == self.params["QUBIC"]["nsub_out"]: + H = self.joint_tod.qubic.H + else: + H = None + + self.joint = JointAcquisitionFrequencyMapMaking( + self.dict_out, + self.params["QUBIC"]["instrument"], + self.params["QUBIC"]["nrec"], + self.params["QUBIC"]["nsub_out"], + H=H, + ) + self.planck_acquisition143 = PlanckAcquisition(143, self.joint.qubic.scene) self.planck_acquisition217 = PlanckAcquisition(217, self.joint.qubic.scene) - self.nus_Q = self._get_averaged_nus() - - ### Joint acquisition for TOD making - if self.params['QUBIC']['nsub_in'] == self.params['QUBIC']['nsub_out']: - self.joint_tod = JointAcquisitionFrequencyMapMaking(self.dict_in, self.params['QUBIC']['instrument'], - self.params['QUBIC']['nsub_in'], - self.params['QUBIC']['nsub_in'], - H=self.joint.qubic.H) - else: - self.joint_tod = JointAcquisitionFrequencyMapMaking(self.dict_in, self.params['QUBIC']['instrument'], - self.params['QUBIC']['nsub_in'], - self.params['QUBIC']['nsub_in'], - H=None) + self.nus_Q = self.get_averaged_nus() ### Coverage map self.coverage = self.joint.qubic.subacqs[0].get_coverage() @@ -150,13 +158,16 @@ def __init__(self, comm, file, params): ### Initial maps self.m_nu_in = self.get_input_map() - + ### Define reconstructed and TOD operator self.get_H() ### Inverse noise covariance matrix - self.invN = self.joint.get_invntt_operator(mask=self.mask) - + if self.params['PLANCK']['external_data']: + self.invN = self.joint.get_invntt_operator(mask=self.mask) + else: + self.invN = self.joint.qubic.get_invntt_operator() + ### Noises seed_noise_planck = self.get_random_value() @@ -217,41 +228,6 @@ def get_components_fgb(self): dict_comps += [Synchrotron(nu0=150, beta_pl=-3)] return dict_comps - - def _get_averaged_nus(self): - - """ - - Method to average QUBIC frequencies. - - """ - - nus_eff = [] - for i in range(self.params['QUBIC']['nrec']): - nus_eff += [np.mean(self.joint.qubic.allnus[i*self.fsub_out:(i+1)*self.fsub_out])] - - return np.array(nus_eff) - - def get_random_value(self): - """Random value - - Method to build a random seed. - - Returns - ------- - seed: int - Random seed. - - """ - - np.random.seed(None) - if self.rank == 0: - seed = np.random.randint(10000000) - else: - seed = None - - seed = self.comm.bcast(seed, root=0) - return seed def get_H(self): """Acquisition operator @@ -269,11 +245,15 @@ def get_H(self): self.H_in_qubic = self.joint_tod.qubic.get_operator(fwhm=self.fwhm_in) ### Pointing matrix for reconstruction - self.H_out_all_pix = self.joint.get_operator(fwhm=self.fwhm_out) - self.H_out = self.joint.get_operator( - fwhm=self.fwhm_out, seenpix=self.seenpix - ) # , external_data=self.params['PLANCK']['external_data']) - + if self.params['PLANCK']['external_data']: + self.H_out_all_pix = self.joint.get_operator(fwhm=self.fwhm_out) + self.H_out = self.joint.get_operator( + fwhm=self.fwhm_out, seenpix=self.seenpix + ) # , external_data=self.params['PLANCK']['external_data']) + else: + self.H_out = self.joint.qubic.get_operator(fwhm=self.fwhm_out) + self.H_out_all_pix = self.joint.qubic.get_operator(fwhm=self.fwhm_out) + def get_averaged_nus(self): """Average frequency @@ -286,11 +266,13 @@ def get_averaged_nus(self): """ - nus_eff = [] - for i in range(self.params['QUBIC']['nrec']): - nus_eff += [np.mean(self.joint.qubic.allnus[i*self.fsub_out:(i+1)*self.fsub_out])] + nus_ave = [] + for i in range(self.params["QUBIC"]["nrec"]): + nus_ave += [ + np.mean(self.joint.qubic.allnus[i * self.fsub_out : (i + 1) * self.fsub_out]) + ] - return np.array(nus_eff) + return np.array(nus_ave) def get_sky_config(self): """Sky configuration. @@ -324,7 +306,7 @@ def get_sky_config(self): seed = self.comm.bcast(seed, root=0) else: seed = self.params["CMB"]["seed"] - print(f"Seed of the CMB is {seed} for rank {self.rank}") + dict_sky["cmb"] = seed for j in self.params["Foregrounds"]: @@ -397,6 +379,27 @@ def get_dict(self, key="in"): dict_qubic[str(i)] = args[i] return dict_qubic + + def get_random_value(self): + """Random value + + Method to build a random seed. + + Returns + ------- + seed: int + Random seed. + + """ + + np.random.seed(None) + if self.rank == 0: + seed = np.random.randint(10000000) + else: + seed = None + + seed = self.comm.bcast(seed, root=0) + return seed def _get_scalar_acquisition_operator(self): """ @@ -416,7 +419,7 @@ def _get_scalar_acquisition_operator(self): acquisition_operators[freq](vector_ones) ) return scalar_acquisition_operators - + def get_convolution(self): """QUBIC resolutions. @@ -434,9 +437,10 @@ def get_convolution(self): """ ### Define FWHMs - fwhm_in = None - fwhm_out = None - + fwhm_in = np.zeros(self.params["QUBIC"]["nsub_in"]) + fwhm_out = np.zeros(self.params["QUBIC"]["nsub_out"]) + fwhm_rec = np.zeros(self.params["QUBIC"]["nrec"]) + ### FWHMs during map-making if self.params["QUBIC"]["convolution_in"]: fwhm_in = self.joint_tod.qubic.allfwhm.copy() @@ -491,13 +495,13 @@ def get_convolution(self): weight_factor = f_cmb.eval(self.joint.qubic.allnus) fun = lambda nu: np.abs(fraction - f_cmb.eval(nu)) - ### Compute expected resolutions and frequencies when not adding convolutions during reconstruction + ### Compute expected resolutions and frequencies when not adding convolutions during reconstruction ### See FMM annexe B to understand the computations for irec in range(self.params["QUBIC"]["nrec"]): numerator_fwhm, denominator_fwhm = 0, 0 numerator_nus, denominator_nus = 0, 0 - for jsub in range(irec * self.fsub_in, (irec + 1) * self.fsub_in): + for jsub in range(irec * self.fsub_out, (irec + 1) * self.fsub_out): # Compute the expected reconstructed resolution for sub-acquisition numerator_fwhm += ( scalar_acquisition_operators[jsub] @@ -547,7 +551,7 @@ def get_input_map(self): m_nu_in = np.zeros( (self.params["QUBIC"]["nrec"], 12 * self.params["SKY"]["nside"] ** 2, 3) ) - + for i in range(self.params["QUBIC"]["nrec"]): m_nu_in[i] = np.mean( self.external_timeline.m_nu[i * self.fsub_out : (i + 1) * self.fsub_out], axis=0 @@ -555,16 +559,11 @@ def get_input_map(self): return m_nu_in - def get_tod(self, noise=False): + def get_tod(self): r"""Simulated TOD. Method that compute observed TODs with :math:`\vec{TOD} = H \cdot \vec{s} + \vec{n}`, with H the QUBIC operator, :math:`\vec{s}` the sky signal and :math:`\vec{n}` the instrumental noise`. - Parameters - ---------- - noise : bool, optional - True if you want to simulate noise in your data, False otherwise, by default False. - Returns ------- TOD: array_like @@ -572,10 +571,6 @@ def get_tod(self, noise=False): """ - if noise: - factor = 0 - else: - factor = 1 if self.params["QUBIC"]["instrument"] == "UWB": if self.params["QUBIC"]["nrec"] != 1: TOD_PLANCK = np.zeros( @@ -589,7 +584,7 @@ def get_tod(self, noise=False): if self.params["QUBIC"]["convolution_in"]: C = HealpixConvolutionGaussianOperator( fwhm=np.min( - self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] + self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in], lmax = 2 * self.params["SKY"]["nside"], ) ) @@ -597,7 +592,7 @@ def get_tod(self, noise=False): C = HealpixConvolutionGaussianOperator(fwhm=0) TOD_PLANCK[irec] = C( - factor * self.external_timeline.maps[irec] + self.noise143 + self.external_timeline.maps[irec] + self.noise143 ) for irec in range( @@ -606,7 +601,7 @@ def get_tod(self, noise=False): if self.params["QUBIC"]["convolution_in"]: C = HealpixConvolutionGaussianOperator( fwhm=np.min( - self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] + self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in], lmax = 2 * self.params["SKY"]["nside"] ) ) @@ -614,7 +609,7 @@ def get_tod(self, noise=False): C = HealpixConvolutionGaussianOperator(fwhm=0) TOD_PLANCK[irec] = C( - factor * self.external_timeline.maps[irec] + self.noise217 + self.external_timeline.maps[irec] + self.noise217 ) else: TOD_PLANCK = np.zeros( @@ -626,20 +621,20 @@ def get_tod(self, noise=False): ) if self.params["QUBIC"]["convolution_in"]: - C = HealpixConvolutionGaussianOperator(fwhm=self.fwhm_in[-1]) + C = HealpixConvolutionGaussianOperator(fwhm=self.fwhm_in[-1], lmax = 2 * self.params["SKY"]["nside"]) else: C = HealpixConvolutionGaussianOperator(fwhm=0) TOD_PLANCK[0] = C( - factor * self.external_timeline.maps[0] + self.noise143 + self.external_timeline.maps[0] + self.noise143 ) TOD_PLANCK[1] = C( - factor * self.external_timeline.maps[0] + self.noise217 + self.external_timeline.maps[0] + self.noise217 ) TOD_PLANCK = TOD_PLANCK.ravel() TOD_QUBIC = ( - self.H_in_qubic(factor * self.external_timeline.m_nu).ravel() + self.H_in_qubic(self.external_timeline.m_nu).ravel() + self.noiseq ) if self.params["PLANCK"]["external_data"]: @@ -649,9 +644,9 @@ def get_tod(self, noise=False): else: - sh_q = self.joint_tod.qubic.ndets * self.joint_tod.qubic.nsamples + sh_q = self.joint.qubic.ndets * self.joint.qubic.nsamples TOD_QUBIC = ( - self.H_in_qubic(factor * self.external_timeline.m_nu).ravel() + self.H_in_qubic(self.external_timeline.m_nu).ravel() + self.noiseq ) if self.params["PLANCK"]["external_data"] == False: @@ -673,7 +668,7 @@ def get_tod(self, noise=False): C = HealpixConvolutionGaussianOperator( fwhm=np.min( self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] - ) + ), lmax = 2 * self.params["SKY"]["nside"] ) else: @@ -682,7 +677,7 @@ def get_tod(self, noise=False): TOD = np.r_[ TOD, C( - factor * self.external_timeline.maps[irec] + self.noise143 + self.external_timeline.maps[irec] + self.noise143 ).ravel(), ] @@ -694,7 +689,7 @@ def get_tod(self, noise=False): C = HealpixConvolutionGaussianOperator( fwhm=np.min( self.fwhm_in[irec * self.fsub_in : (irec + 1) * self.fsub_in] - ) + ), lmax = 2 * self.params["SKY"]["nside"] ) else: @@ -703,22 +698,12 @@ def get_tod(self, noise=False): TOD = np.r_[ TOD, C( - factor * self.external_timeline.maps[irec] + self.noise217 + self.external_timeline.maps[irec] + self.noise217 ).ravel(), ] return TOD - def _barrier(self): - """ - Method to introduce comm._Barrier() function if MPI communicator is detected. - - """ - if self.comm is None: - pass - else: - self.comm.Barrier() - def _print_message(self, message): """ Method to print message only on rank 0 if MPI communicator is detected. It display simple message if not. @@ -773,8 +758,10 @@ def get_preconditioner(self): ) conditioner[conditioner == np.inf] = 1 - - M = DiagonalOperator(conditioner[:, self.seenpix, :]) + if self.params['PLANCK']['external_data']: + M = DiagonalOperator(conditioner[:, self.seenpix, :]) + else: + M = DiagonalOperator(conditioner) else: M = None return M @@ -805,12 +792,13 @@ def pcg(self, d, x0, seenpix): ### Update components when pixels outside the patch are fixed (assumed to be 0) A = self.H_out.T * self.invN * self.H_out - x_planck = self.m_nu_in * (1 - seenpix[None, :, None]) - print(self.H_out_all_pix.shapein) - print(x_planck.shape) - print(d.shape) - print(self.invN) - b = self.H_out.T * self.invN * (d - self.H_out_all_pix(x_planck)) + if self.params['PLANCK']['external_data']: + x_planck = self.m_nu_in * (1 - seenpix[None, :, None]) + b = self.H_out.T * self.invN * (d - self.H_out_all_pix(x_planck)) + else: + R = ReshapeOperator(shapein=d.shape, shapeout=self.invN.shapein) + b = self.H_out.T * self.invN(R) * d + ### Preconditionning M = self.get_preconditioner() @@ -818,7 +806,13 @@ def pcg(self, d, x0, seenpix): gif_folder = self.plot_folder + f"{self.job_id}/iter/" else: gif_folder = None - + + true_maps = self.m_nu_in.copy() + + for irec in range(self.params["QUBIC"]["nrec"]): + C = HealpixConvolutionGaussianOperator(fwhm=self.fwhm_rec[irec], lmax = 2 * self.joint.qubic.scene.nside) + true_maps[irec] = C(self.m_nu_in[irec]) + ### PCG solution_qubic_planck = pcg( A=A, @@ -836,7 +830,8 @@ def pcg(self, d, x0, seenpix): center=self.center, reso=self.params["PCG"]["resolution_plot"], fwhm_plot=self.params["PCG"]["fwhm_plot"], - input=self.m_nu_in, + input=true_maps, + is_planck=self.params['PLANCK']['external_data'], ) if self.params["PCG"]["gif"]: @@ -850,7 +845,10 @@ def pcg(self, d, x0, seenpix): ) solution = np.ones(self.m_nu_in.shape) * hp.UNSEEN - solution[:, seenpix, :] = solution_qubic_planck["x"]["x"].copy() + if self.params['PLANCK']['external_data']: + solution[:, seenpix, :] = solution_qubic_planck["x"]["x"].copy() + else: + solution[:, seenpix, :] = solution_qubic_planck["x"]["x"][:, seenpix, :].copy() return solution @@ -862,6 +860,16 @@ def _save_data(self, name, d): with open(name, "wb") as handle: pickle.dump(d, handle, protocol=pickle.HIGHEST_PROTOCOL) + + def _barrier(self): + """ + Method to introduce comm._Barrier() function if MPI communicator is detected. + + """ + if self.comm is None: + pass + else: + self.comm.Barrier() def run(self): """Run the FMM Pipeline. @@ -873,21 +881,27 @@ def run(self): self._print_message("\n=========== Map-Making ===========\n") ### Get simulated data - self.TOD = self.get_tod(noise=False) - self.n = self.get_tod(noise=True) + self.TOD = self.get_tod() ### Wait for all processes self._barrier() - ### Build starting point for PCG - starting_maps = self.m_nu_in.copy() - if self.params['QUBIC']['convolution_in']: - for i in range(self.params['QUBIC']['nrec']): - C_rec = HealpixConvolutionGaussianOperator(fwhm=self.fwhm_rec[i]) - starting_maps[i] = C_rec(self.m_nu_in[i]) + #maps_in = self.m_nu_in.copy() + #if self.params["QUBIC"]["convolution_in"]: + # for i in range(self.params["QUBIC"]["nrec"]): + # C = HealpixConvolutionGaussianOperator(fwhm=self.fwhm_rec[i]) + # maps_in[i] = C(self.m_nu_in[i]) + #x0 = self.m_nu_in.copy() + #x0[:, self.seenpix, :] = 0 + + if self.params['PLANCK']['external_data']: + starting_point = np.zeros(self.m_nu_in[:, self.seenpix, :].shape) + else: + starting_point = np.zeros(self.m_nu_in.shape) + ### Solve map-making equation - self.s_hat = self.pcg(self.TOD, x0=starting_maps, seenpix=self.seenpix) + self.s_hat = self.pcg(self.TOD, x0=starting_point, seenpix=self.seenpix) ### Wait for all processes self._barrier() @@ -907,7 +921,8 @@ def run(self): self.external_maps_noise = self.externaldata_noise.maps.copy() self.external_maps_noise[:, ~self.seenpix, :] = 0 - + + self.nus_rec = self.nus_Q.copy() if len(self.externaldata.external_nus) != 0: fwhm_ext = self.externaldata.fwhm_ext.copy() self.s_hat = np.concatenate((self.s_hat, self.external_maps), axis=0) @@ -928,14 +943,7 @@ def run(self): filename=self.plot_folder + f"/all_maps.png", figsize=(10, 5), ) - # self.plots.plot_FMM(self.m_nu_in, self.s_hat, self.center, self.seenpix, self.nus_Q, job_id=self.job_id, istk=0, nsig=3, name='signal') - # self.plots.plot_FMM(self.m_nu_in, self.s_hat, self.center, self.seenpix, self.nus_Q, job_id=self.job_id, istk=1, nsig=3, name='signal') - # self.plots.plot_FMM(self.m_nu_in, self.s_hat, self.center, self.seenpix, self.nus_Q, job_id=self.job_id, istk=2, nsig=3, name='signal') - - # self.plots.plot_FMM(self.m_nu_in*0, self.s_hat_noise, self.center, self.seenpix, self.nus_Q, job_id=self.job_id, istk=0, nsig=3, name='noise') - # self.plots.plot_FMM(self.m_nu_in*0, self.s_hat_noise, self.center, self.seenpix, self.nus_Q, job_id=self.job_id, istk=1, nsig=3, name='noise') - # self.plots.plot_FMM(self.m_nu_in*0, self.s_hat_noise, self.center, self.seenpix, self.nus_Q, job_id=self.job_id, istk=2, nsig=3, name='noise') - + mapmaking_time = time.time() - self.mapmaking_time_0 if self.comm is None: print(f"Map-making done in {mapmaking_time:.3f} s") @@ -944,6 +952,7 @@ def run(self): print(f"Map-making done in {mapmaking_time:.3f} s") dict_solution = { + "maps_in":self.m_nu_in, "maps": self.s_hat, "maps_noise": self.s_hat_noise, "nus": self.nus_rec, @@ -958,7 +967,9 @@ def run(self): } self._save_data(self.file, dict_solution) - self._barrier() + + ### Wait for all processors + _barrier() class PipelineEnd2End: @@ -1029,4 +1040,4 @@ def main(self, specific_file=None): "parameters": self.params, } - self._save_data(self.file_spectrum, dict_solution) + self._save_data(self.file_spectrum, dict_solution) \ No newline at end of file diff --git a/src/lib/Qcg.py b/src/lib/Qcg.py index 42d8144..b101e2b 100644 --- a/src/lib/Qcg.py +++ b/src/lib/Qcg.py @@ -43,6 +43,7 @@ def __init__( input=None, fwhm=0, iter_init=0, + is_planck=False, ): """ Parameters @@ -98,6 +99,7 @@ def __init__( self.fwhm_plot = fwhm_plot self.input = input self.fwhm = fwhm_plot + self.is_planck = is_planck if self.gif is not None: if not os.path.isdir(self.gif): @@ -180,35 +182,35 @@ def iteration(self): self.A(self.d, self.q) alpha = self.delta / self.dot(self.d, self.q) self.x += alpha * self.d - print(self.x.shape) - print(self.input[:, self.seenpix, :].shape) - _r = self.x - self.input[:, self.seenpix, :] + + map_i = self.x.copy() + if self.is_planck: + map_i = np.ones(self.input.shape) * hp.UNSEEN + map_i[:, self.seenpix, :] = self.x.copy() + + _r = map_i[:, self.seenpix, :] - self.input[:, self.seenpix, :] self.rms = np.std(_r, axis=1) + if self.gif is not None: if self.comm.Get_rank() == 0: - mymap = np.zeros(self.input.shape) nsig = 2 min, max = -nsig * np.std( self.input[0, self.seenpix], axis=0 ), nsig * np.std(self.input[0, self.seenpix], axis=0) - mymap[:, self.seenpix, :] = self.x.copy() - mymap[:, ~self.seenpix_plot, :] = hp.UNSEEN - _plot_reconstructed_maps( - mymap, + map_i, self.input, self.gif + f"iter_{self.niterations+self.iter_init}.png", self.center, reso=self.reso, - figsize=(12, 2.7*mymap.shape[0]), + figsize=(12, 2.7*map_i.shape[0]), min=min, max=max, fwhm=self.fwhm, iter=self.niterations, ) - self.r -= alpha * self.q self.error = np.sqrt(self.norm(self.r) / self.b_norm) self.convergence = np.append(self.convergence, self.error) @@ -252,6 +254,7 @@ def pcg( input=None, fwhm=0, iter_init=0, + is_planck=False, ): """ output = pcg(A, b, [x0, tol, maxiter, M, disp, callback, @@ -321,6 +324,7 @@ def pcg( input=input, fwhm=fwhm, iter_init=iter_init, + is_planck=is_planck, ) try: output = algo.run() @@ -353,4 +357,4 @@ def _dot(x, y, comm): d = np.array(np.dot(x.ravel(), y.ravel())) if comm is not None: comm.Allreduce(MPI.IN_PLACE, d) - return d + return d \ No newline at end of file diff --git a/src/lib/Qmap_plotter.py b/src/lib/Qmap_plotter.py index e43f228..a445326 100644 --- a/src/lib/Qmap_plotter.py +++ b/src/lib/Qmap_plotter.py @@ -809,4 +809,4 @@ def plot_rms_iteration(self, rms, figsize=(8, 6), ki=0): if ki > 0: os.remove(f"CMM/jobs/{self.job_id}/rms_iter{ki}.png") - plt.close() + plt.close() \ No newline at end of file