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..6a00fa2 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: 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 bandpass_correction: True # Correction to keep the edges of the integration band constant NOISE: @@ -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 @@ -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,7 +60,7 @@ 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 gif: True diff --git a/src/FMM/pipeline.py b/src/FMM/pipeline.py index dc218a4..728f0f5 100644 --- a/src/FMM/pipeline.py +++ b/src/FMM/pipeline.py @@ -22,6 +22,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 @@ -35,7 +36,7 @@ def save_pkl(name, d): """ with open(name, "wb") as handle: pickle.dump(d, handle, protocol=pickle.HIGHEST_PROTOCOL) - +''' __all__ = ["PipelineFrequencyMapMaking", "PipelineEnd2End"] @@ -66,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/" @@ -147,7 +149,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"], @@ -160,8 +162,11 @@ def __init__(self, comm, file, params): 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() @@ -223,27 +228,6 @@ def get_components_fgb(self): return dict_comps - 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 @@ -260,11 +244,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 @@ -280,7 +268,7 @@ 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]) + np.mean(self.joint.qubic.allnus[i * self.fsub_out : (i + 1) * self.fsub_out]) ] return np.array(nus_ave) @@ -317,7 +305,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"]: @@ -390,6 +378,47 @@ 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): + """ + Function that will compute "scalar acquisition operatord" by applying the acquisition operators to a vector full of ones. + These scalar operators will be used to compute the resolutions in the case where we do not add convolutions during reconstruction. + """ + ### Import the acquisition operators + acquisition_operators = self.joint.qubic.H + + ### Create the vector full of ones which will be used to compute the scalar operators + vector_ones = np.ones(acquisition_operators[0].shapein) + + ### Apply each sub_operator on the vector + scalar_acquisition_operators = np.empty(len(self.joint.qubic.allnus)) + for freq in range(len(self.joint.qubic.allnus)): + scalar_acquisition_operators[freq] = np.mean( + acquisition_operators[freq](vector_ones) + ) + return scalar_acquisition_operators def get_convolution(self): """QUBIC resolutions. @@ -408,9 +437,9 @@ 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"]: @@ -455,7 +484,41 @@ def get_convolution(self): and self.params["QUBIC"]["convolution_out"] is False ): fwhm_rec = np.array([]) + scalar_acquisition_operators = self._get_scalar_acquisition_operator() + + if self.params["Foregrounds"]["Dust"]: + f_dust = c.ModifiedBlackBody(nu0=353, beta_d=1.54) + weight_factor = f_dust.eval(self.joint.qubic.allnus) + fun = lambda nu: np.abs(fraction - f_dust.eval(nu)) + else: + f_cmb = c.CMB() + 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 + ### 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_out, (irec + 1) * self.fsub_out): + # Compute the expected reconstructed resolution for sub-acquisition + numerator_fwhm += ( + scalar_acquisition_operators[jsub] + * weight_factor[jsub] + * fwhm_in[jsub] + ) + denominator_fwhm += ( + scalar_acquisition_operators[jsub] * weight_factor[jsub] + ) + + # Compute the expected reconstructed frequencies for sub_acquisition + numerator_nus += ( + scalar_acquisition_operators[jsub] * weight_factor[jsub] + ) + denominator_nus += scalar_acquisition_operators[jsub] + + # Compute the expected resolution fwhm_rec = np.append( fwhm_rec, np.mean( @@ -495,21 +558,16 @@ 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 - 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 @@ -517,10 +575,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( @@ -534,7 +588,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], lmax = 2 * self.params["SKY"]["nside"], ) ) @@ -542,7 +596,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( @@ -551,7 +605,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], lmax = 2 * self.params["SKY"]["nside"] ) ) @@ -559,7 +613,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( @@ -571,20 +625,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"]: @@ -596,7 +650,7 @@ def get_tod(self, noise=False): 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: @@ -617,8 +671,8 @@ 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] + ), lmax = 2 * self.params["SKY"]["nside"] ) else: @@ -627,7 +681,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(), ] @@ -638,8 +692,8 @@ 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] + ), lmax = 2 * self.params["SKY"]["nside"] ) else: @@ -648,22 +702,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. @@ -710,16 +754,18 @@ 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) ) 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 @@ -750,8 +796,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]) - 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() @@ -759,7 +810,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, @@ -777,7 +834,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"]: @@ -791,7 +849,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 @@ -803,6 +864,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. @@ -814,21 +885,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() @@ -848,7 +925,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) @@ -869,14 +947,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") @@ -885,6 +956,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, @@ -899,7 +971,9 @@ def run(self): } self._save_data(self.file, dict_solution) - self._barrier() + + ### Wait for all processors + _barrier() class PipelineEnd2End: @@ -970,4 +1044,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 75812f0..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,33 +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 - _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) @@ -250,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, @@ -319,6 +324,7 @@ def pcg( input=input, fwhm=fwhm, iter_init=iter_init, + is_planck=is_planck, ) try: output = algo.run() @@ -351,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 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: