From e70f6f07da89e3743e5adcf77d88fe85a90fd532 Mon Sep 17 00:00:00 2001 From: Reijo Keskitalo Date: Fri, 22 Apr 2022 12:36:12 -0700 Subject: [PATCH 1/3] Add support to HWP in beam convolution --- cmake.sh | 1 + src/toast/ops/__init__.py | 2 +- src/toast/ops/conviqt.py | 137 +++++++- src/toast/tests/_helpers.py | 81 +++-- src/toast/tests/ops_sim_tod_conviqt.py | 462 ++++++++++++++++++++----- 5 files changed, 554 insertions(+), 129 deletions(-) create mode 100755 cmake.sh diff --git a/cmake.sh b/cmake.sh new file mode 100755 index 000000000..450569532 --- /dev/null +++ b/cmake.sh @@ -0,0 +1 @@ +cmake .. -DCMAKE_INSTALL_PREFIX=$PREFIX diff --git a/src/toast/ops/__init__.py b/src/toast/ops/__init__.py index b5a456543..d56c7d8f5 100644 --- a/src/toast/ops/__init__.py +++ b/src/toast/ops/__init__.py @@ -61,7 +61,7 @@ from .madam import Madam, madam_params_from_mapmaker -from .conviqt import SimConviqt, SimWeightedConviqt +from .conviqt import SimConviqt, SimWeightedConviqt, SimTEBConviqt from .sim_gaindrifts import GainDrifter from .sim_crosstalk import CrossTalk, MitigateCrossTalk diff --git a/src/toast/ops/conviqt.py b/src/toast/ops/conviqt.py index eec783f49..fe3c39152 100644 --- a/src/toast/ops/conviqt.py +++ b/src/toast/ops/conviqt.py @@ -266,6 +266,8 @@ def _exec(self, data, detectors=None, **kwargs): detector = self.get_detector(det) theta, phi, psi, psi_pol = self.get_pointing(data, det, verbose) + np.savez(f'{self.det_data}_{det}.npz', psi=psi, psi_pol=psi_pol, + theta=theta, phi=phi) pnt = self.get_buffer(theta, phi, psi, det, verbose) del theta, phi, psi @@ -610,13 +612,8 @@ def _exec(self, data, detectors=None, **kwargs): beam_file = self.beam_file_dict[det] else: beam_file = self.beam_file.format(detector=det, mc=self.mc) - beam_file_i00 = beam_file.replace(".fits", "_I000.fits") - beam_file_0i0 = beam_file.replace(".fits", "_0I00.fits") - beam_file_00i = beam_file.replace(".fits", "_00I0.fits") - beamI00 = self.get_beam(beam_file_i00, det, verbose) - beam0I0 = self.get_beam(beam_file_0i0, det, verbose) - beam00I = self.get_beam(beam_file_00i, det, verbose) + beamI00, beam0I0, beam00I = self.get_beam(beam_file, det, verbose) detector = self.get_detector(det) @@ -656,3 +653,131 @@ def _exec(self, data, detectors=None, **kwargs): timer.report_clear(f"conviqt process detector {det}") return + + def get_beam(self, beamfile, det, verbose): + timer = Timer() + timer.start() + beam_file_i00 = beamfile.replace(".fits", "_I000.fits") + beam_file_0i0 = beamfile.replace(".fits", "_0I00.fits") + beam_file_00i = beamfile.replace(".fits", "_00I0.fits") + beami00 = conviqt.Beam( + self.lmax, self.beammmax, self.pol, beam_file_i00, self.comm + ) + beam0i0 = conviqt.Beam( + self.lmax, self.beammmax, self.pol, beam_file_0i0, self.comm + ) + beam00i = conviqt.Beam( + self.lmax, self.beammmax, self.pol, beam_file_00i, self.comm + ) + + if verbose: + timer.report_clear(f"initialize beam for detector {det}") + return beami00, beam0i0, beam00i + + +class SimTEBConviqt(SimConviqt): + """ + Operator that uses libconviqt to generate beam-convolved timestreams. + This operator should be used in presence of a spinning HWP which makes the beam time-dependent, + constantly mapping the co- and cross-polar responses on to each other. + In the parent class OpSimConviqt we assume the beam to be static. + + + The convolution is performed by coupling each IQU component of the signal propertly as: + :math:`skyT_lm * beamT_lm, skyE_lm * Re{P}, skyB_lm * Im{P}`. + + For extra details please refer to [this note ](https://giuspugl.github.io/reports/Notes_TEB_convolution.html) + """ + + @function_timer + def _exec(self, data, detectors=None, **kwargs): + if not self.available: + raise RuntimeError("libconviqt is not available") + + if self.comm is None: + raise RuntimeError("libconviqt requires MPI") + + if self.detector_pointing is None: + raise RuntimeError("detector_pointing cannot be None.") + + log = Logger.get() + + timer = Timer() + timer.start() + + # Expand detector pointing + self.detector_pointing.apply(data, detectors=detectors) + + all_detectors = self._get_all_detectors(data, detectors) + + for det in all_detectors: + verbose = self.comm.rank == 0 and self.verbosity > 0 + + # Expand detector pointing + self.detector_pointing.apply(data, detectors=[det]) + + if det in self.sky_file_dict: + sky_file = self.sky_file_dict[det] + else: + sky_file = self.sky_file.format(detector=det, mc=self.mc) + skyT = self.get_sky(sky_file.replace(".fits","_T.fits"), det, verbose) + skyE = self.get_sky(sky_file.replace(".fits","_E.fits"), det, verbose) + skyB = self.get_sky(sky_file.replace(".fits","_B.fits"), det, verbose) + + if det in self.beam_file_dict: + beam_file = self.beam_file_dict[det] + else: + beam_file = self.beam_file.format(detector=det, mc=self.mc) + + beam_T, beam_P= self.get_beam(beam_file, det, verbose) + + detector = self.get_detector(det) + + + theta, phi, psi, psi_pol = self.get_pointing(data, det, verbose) + # T-convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + + convolved_data = self.convolve(skyT, beam_T, detector, pnt, det, verbose) + + del (pnt,) + # E-convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.cos(2 * psi_pol) * self.convolve( + skyE, beam_P, detector, pnt, det, verbose + ) + del (pnt,) + # B-convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.sin(2 * psi_pol) * self.convolve( + skyB, beam_P, detector, pnt, det, verbose + ) + del theta, phi, psi + + self.calibrate_signal( + data, + det, + beam_T, + convolved_data, + verbose, + ) + self.save(data, det, convolved_data, verbose) + + del pnt, detector, beam_T, beam_P , skyT, skyE, skyB + + if verbose: + timer.report_clear(f"conviqt process detector {det}") + + return + + def get_beam(self, beamfile, det, verbose): + timer = Timer() + timer.start() + beam_file_T = beamfile.replace(".fits", "_T.fits") + beam_file_P = beamfile.replace(".fits", "_P.fits") + beamT = conviqt.Beam(lmax = self.lmax, mmax = self.beammmax, pol = self.pol , beamfile=beam_file_T, comm= self.comm) + beamP = conviqt.Beam(lmax = self.lmax, mmax = self.beammmax, pol = self.pol , beamfile=beam_file_P, comm= self.comm) + + if verbose: + timer.report_clear(f"initialize beam for detector {det}") + return beamT, beamP diff --git a/src/toast/tests/_helpers.py b/src/toast/tests/_helpers.py index 0c28893c3..bceed1557 100644 --- a/src/toast/tests/_helpers.py +++ b/src/toast/tests/_helpers.py @@ -521,52 +521,71 @@ def create_fake_beam_alm( mmax=10, fwhm_x=10 * u.degree, fwhm_y=10 * u.degree, - pol=True, + pol = True, separate_IQU=False, + separate_TP=False , + detB_beam=False, + normalize_beam=False, ): # pick an nside >= lmax to be sure that the a_lm will be fairly accurate nside = 2 while nside < lmax: - nside *= 2 npix = 12 * nside**2 pix = np.arange(npix) - vec = hp.pix2vec(nside, pix, nest=False) - theta, phi = hp.vec2dir(vec) - x = theta * np.cos(phi) - y = theta * np.sin(phi) - sigma_x = fwhm_x.to_value(u.radian) / np.sqrt(8 * np.log(2)) + x, y, z = hp.pix2vec(nside, pix, nest=False) + sigma_z = fwhm_x.to_value(u.radian) / np.sqrt(8 * np.log(2)) sigma_y = fwhm_y.to_value(u.radian) / np.sqrt(8 * np.log(2)) - beam_map = np.exp(-0.5 * (x**2 / sigma_x**2 + y**2 / sigma_y**2)) - empty = np.zeros_like(beam_map) - if pol and separate_IQU: - beam_map_I = np.vstack([beam_map, empty, empty]) - beam_map_Q = np.vstack([empty, beam_map, empty]) - beam_map_U = np.vstack([empty, empty, beam_map]) + beam = np.exp(-((z ** 2 / 2 / sigma_z ** 2 + y ** 2 / 2 / sigma_y ** 2))) + beam[x < 0] = 0 + beam_map = np.zeros([3, npix]) + beam_map[0] = beam + beam_map[1] = beam + bl, blm = hp.anafast(beam_map, lmax=lmax, iter=0, alm=True, pol=pol) + hp.rotate_alm(blm, psi=0, theta=-np.pi / 2, phi=0) + + if detB_beam: + # we make sure that the two detectors within the same pair encode + # two beams with the flipped sign in Q U beams + beam_map = hp.alm2map(blm, nside=nside, lmax=lmax, mmax=mmax) + beam_map[1:] *= -1 + blm = hp.map2alm(beam_map, lmax=lmax, mmax=mmax) + + if normalize_beam: + # We make sure that the simulated beams are normalized in the test + # for the normalization we follow the convention adopted in Conviqt, + # i.e. the monopole term in the map is left unchanged + idx = hp.Alm.getidx(lmax=lmax, l=0, m=0) + norm = 2 * np.pi * blm[0, idx].real - try: - a_lm = [ - hp.map2alm(beam_map_I, lmax=lmax, mmax=mmax, verbose=False), - hp.map2alm(beam_map_Q, lmax=lmax, mmax=mmax, verbose=False), - hp.map2alm(beam_map_U, lmax=lmax, mmax=mmax, verbose=False), - ] - except TypeError: - # older healpy which does not have verbose keyword - a_lm = [ - hp.map2alm(beam_map_I, lmax=lmax, mmax=mmax), - hp.map2alm(beam_map_Q, lmax=lmax, mmax=mmax), - hp.map2alm(beam_map_U, lmax=lmax, mmax=mmax), - ] else: - if pol: - beam_map = np.vstack([beam_map, beam_map, empty]) + norm=1. + + blm /= norm + if separate_IQU: + empty = np.zeros_like(beam_map[0]) + beam_map_I = np.vstack([beam_map[0], empty, empty]) + beam_map_Q = np.vstack([empty, beam_map[1], empty]) + beam_map_U = np.vstack([empty, empty, beam_map[1]]) try: - a_lm = hp.map2alm(beam_map, lmax=lmax, mmax=mmax, verbose=False) + blmi00 = hp.map2alm(beam_map_I, lmax=lmax, mmax=mmax, verbose=False, pol=True)/norm + blm0i0 = hp.map2alm(beam_map_Q, lmax=lmax, mmax=mmax, verbose=False, pol=True )/norm + blm00i= hp.map2alm(beam_map_U, lmax=lmax, mmax=mmax, verbose=False, pol=True)/norm except TypeError: - a_lm = hp.map2alm(beam_map, lmax=lmax, mmax=mmax) + # older healpy which does not have verbose keyword + blmi00 = hp.map2alm(beam_map_I, lmax=lmax, mmax=mmax,pol=True)/norm + blm0i0 = hp.map2alm(beam_map_Q, lmax=lmax, mmax=mmax,pol=True)/norm + blm00i = hp.map2alm(beam_map_U, lmax=lmax, mmax=mmax,pol=True)/norm + return [blmi00,blm0i0,blm00i] - return a_lm + elif separate_TP : + blmT =blm[0].copy() + blmP =blm[1 ] + blm[2] + + return [blmT, blmP ] + else: + return blm def fake_flags( diff --git a/src/toast/tests/ops_sim_tod_conviqt.py b/src/toast/tests/ops_sim_tod_conviqt.py index 9121398f2..3c56d8be2 100644 --- a/src/toast/tests/ops_sim_tod_conviqt.py +++ b/src/toast/tests/ops_sim_tod_conviqt.py @@ -33,13 +33,15 @@ class SimConviqtTest(MPITestCase): def setUp(self): + + np.random.seed(777) fixture_name = os.path.splitext(os.path.basename(__file__))[0] self.outdir = create_outdir(self.comm, fixture_name) - self.nside = 64 + self.nside = 32 self.lmax = 128 self.fwhm_sky = 10 * u.degree - self.fwhm_beam = 30 * u.degree + self.fwhm_beam = 15 * u.degree self.mmax = self.lmax self.fname_sky = os.path.join(self.outdir, "sky_alm.fits") self.fname_beam = os.path.join(self.outdir, "beam_alm.fits") @@ -51,15 +53,29 @@ def setUp(self): if self.rank == 0: # Synthetic sky and beam (a_lm expansions) self.slm = create_fake_sky_alm(self.lmax, self.fwhm_sky) - self.slm[1:] = 0 # No polarization - hp.write_alm(self.fname_sky, self.slm, lmax=self.lmax, overwrite=True) + hp.write_alm(self.fname_sky, self.slm, lmax=self.lmax, overwrite=True) + hp.write_alm(self.fname_sky.replace(".fits","_T.fits"), self.slm[0] , lmax=self.lmax, overwrite=True) + hp.write_alm(self.fname_sky.replace(".fits","_E.fits"), self.slm[1] , lmax=self.lmax, overwrite=True) + hp.write_alm(self.fname_sky.replace(".fits","_B.fits"), self.slm[2] , lmax=self.lmax, overwrite=True) + self.blm = create_fake_beam_alm( self.lmax, self.mmax, fwhm_x=self.fwhm_beam, fwhm_y=self.fwhm_beam, + normalize_beam=True, + ) + + self.blm_bottom = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + normalize_beam=True, + detB_beam=True, ) + hp.write_alm( self.fname_beam, self.blm, @@ -67,13 +83,29 @@ def setUp(self): mmax_in=self.mmax, overwrite=True, ) - + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom.fits"), + self.blm_bottom, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) blm_I, blm_Q, blm_U = create_fake_beam_alm( self.lmax, self.mmax, fwhm_x=self.fwhm_beam, fwhm_y=self.fwhm_beam, separate_IQU=True, + normalize_beam=True, + ) + blm_Ibot, blm_Qbot, blm_Ubot = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + separate_IQU=True, + detB_beam=True, + normalize_beam=True, ) hp.write_alm( self.fname_beam.replace(".fits", "_I000.fits"), @@ -96,19 +128,111 @@ def setUp(self): mmax_in=self.mmax, overwrite=True, ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_I000.fits"), + blm_Ibot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_0I00.fits"), + blm_Qbot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_00I0.fits"), + blm_Ubot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + + # we explicitly store 3 separate beams for the T, E and B sky alm. + blm_T, blm_P = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + separate_TP=True, + detB_beam=False, + normalize_beam=True, + ) + blm_Tbot, blm_Pbot = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + separate_TP=True, + detB_beam=True, + normalize_beam=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_T.fits"), + blm_T, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_T.fits"), + blm_Tbot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_P.fits"), + blm_P, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_P.fits"), + blm_Pbot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + # in if self.comm is not None: self.comm.barrier() return - def test_sim(self): + def make_beam_file_dict(self, data): + """ + We make sure that data observed in each A/B detector within a + detector pair will have the right beam associated to it. In particular, + Q and U beams of B detector must have a flipped sign wrt the A one. + """ + + fname2 = self.fname_beam.replace(".fits", "_bottom.fits") + + self.beam_file_dict = {} + for det in data.obs[0].local_detectors: + if det[-1] == "A": + self.beam_file_dict[det] = self.fname_beam + else: + self.beam_file_dict[det] = fname2 + + return + + def test_sim_conviqt(self): if not ops.conviqt.available(): print("libconviqt not available, skipping tests") return # Create a fake scan strategy that hits every pixel once. - data = create_healpix_ring_satellite(self.comm, nside=self.nside) + # data = create_healpix_ring_satellite(self.comm, nside=self.nside) + data = create_satellite_data( + self.comm, obs_time=120 * u.min, pixel_per_process=2 + ) + self.make_beam_file_dict(data) # Generate timestreams @@ -119,12 +243,14 @@ def test_sim(self): comm=self.comm, detector_pointing=detpointing, sky_file=self.fname_sky, - beam_file=self.fname_beam, + beam_file_dict=self.beam_file_dict, dxx=False, det_data=key, - normalize_beam=True, + pol=True, + normalize_beam=False, # beams are already produced normalized fwhm=self.fwhm_sky, ) + sim_conviqt.apply(data) # Bin a map to study @@ -136,7 +262,8 @@ def test_sim(self): ) pixels.apply(data) weights = ops.StokesWeights( - mode="I", + mode="IQU", + hwp_angle=None, detector_pointing=detpointing, ) weights.apply(data) @@ -177,23 +304,14 @@ def test_sim(self): fail = False if self.rank == 0: - set_matplotlib_backend() - - import matplotlib.pyplot as plt - hitsfile = os.path.join(self.outdir, "toast_hits.fits") + hdata = hp.read_map(hitsfile) - hp.mollview(hdata, xsize=1600) - plt.savefig(os.path.join(self.outdir, "toast_hits.png")) - plt.close() + footprint = np.ma.masked_not_equal(hdata, 0.0).mask mapfile = os.path.join(self.outdir, "toast_bin.fits") - mdata = hp.read_map(mapfile) - - hp.mollview(mdata, xsize=1600) - plt.savefig(os.path.join(self.outdir, "toast_bin.png")) - plt.close() + mdata = hp.read_map(mapfile, field=range(3)) deconv = 1 / hp.gauss_beam( self.fwhm_sky.to_value(u.radian), @@ -202,99 +320,263 @@ def test_sim(self): ) smoothed = hp.alm2map( - hp.almxfl(self.slm[0], deconv), + [hp.almxfl(self.slm[ii], deconv) for ii in range(3)], self.nside, lmax=self.lmax, fwhm=self.fwhm_beam.to_value(u.radian), verbose=False, pixwin=False, ) - + smoothed[:, ~footprint] = 0 cl_out = hp.anafast(mdata, lmax=self.lmax) cl_smoothed = hp.anafast(smoothed, lmax=self.lmax) - cl_in = hp.alm2cl(self.slm[0], lmax=self.lmax) - blsq = hp.alm2cl(self.blm[0], lmax=self.lmax, mmax=self.mmax) - blsq /= blsq[1] + np.testing.assert_almost_equal(cl_smoothed[0], cl_out[0], decimal=2) - gauss_blsq = hp.gauss_beam( - self.fwhm_beam.to_value(u.radian), - lmax=self.lmax, - pol=False, - ) + return - sky = hp.alm2map(self.slm[0], self.nside, lmax=self.lmax, verbose=False) - beam = hp.alm2map( - self.blm[0], - self.nside, - lmax=self.lmax, - mmax=self.mmax, - verbose=False, - ) + def test_sim_weighted_conviqt(self): + if not ops.conviqt.available(): + print("libconviqt not available, skipping tests") + return - fig = plt.figure(figsize=[12, 8]) - nrow, ncol = 2, 3 - hp.mollview(sky, title="input sky", sub=[nrow, ncol, 1]) - hp.mollview(mdata, title="output sky", sub=[nrow, ncol, 2]) - hp.mollview(smoothed, title="smoothed sky", sub=[nrow, ncol, 3]) - hp.mollview(beam, title="beam", sub=[nrow, ncol, 4], rot=[0, 90]) - - ell = np.arange(self.lmax + 1) - ax = fig.add_subplot(nrow, ncol, 5) - ax.plot(ell[1:], cl_in[1:], label="input") - ax.plot(ell[1:], cl_smoothed[1:], label="smoothed") - ax.plot(ell[1:], blsq[1:], label="beam") - ax.plot(ell[1:], gauss_blsq[1:], label="gauss beam") - ax.plot(ell[1:], 1 / deconv[1:] ** 2, label="1 / deconv") - ax.plot( - ell[1:], - cl_in[1:] * blsq[1:] * deconv[1:] ** 2, - label="input x beam x deconv", - ) - ax.plot(ell[1:], cl_out[1:], label="output") - ax.legend(loc="best") - ax.set_xscale("log") - ax.set_yscale("log") - ax.set_ylim([1e-20, 1e1]) - - # For some reason, matplotlib hangs with multiple tasks, - # even if only one writes. Uncomment these lines when debugging. - # - # outfile = os.path.join(self.outdir, "cl_comparison.png") - # plt.savefig(outfile) - # plt.close() - - compare = blsq > 1e-5 - ref = cl_in[compare] * blsq[compare] * deconv[compare] ** 2 - norm = np.mean(cl_out[compare] / ref) - if not np.allclose( - norm * ref, - cl_out[compare], - rtol=1e-5, - atol=1e-5, - ): - fail = True + # Create a fake scan strategy that hits every pixel once. + # data = create_healpix_ring_satellite(self.comm, nside=self.nside) + data = create_satellite_data( + self.comm, obs_time=120 * u.min, pixel_per_process=2 + ) + self.make_beam_file_dict(data) - if self.comm is not None: - fail = self.comm.bcast(fail, root=0) + # Generate timestreams + + detpointing = ops.PointingDetectorSimple() + + key1 = "conviqt" + sim_conviqt = ops.SimConviqt( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False, + det_data=key1, + pol=True, + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + sim_conviqt.apply(data) + + key2 = "wconviqt" + + sim_wconviqt = ops.SimWeightedConviqt( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False, + det_data=key2, + pol=True, + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + sim_wconviqt.apply(data) + # Bin a map to study + + pixels = ops.PixelsHealpix( + nside=self.nside, + nest=False, + detector_pointing=detpointing, + ) + pixels.apply(data) + weights = ops.StokesWeights( + mode="IQU", + hwp_angle=None, + detector_pointing=detpointing, + ) + weights.apply(data) - self.assertFalse(fail) + default_model = ops.DefaultNoiseModel() + default_model.apply(data) + + cov_and_hits = ops.CovarianceAndHits( + pixel_dist="pixel_dist", + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + rcond_threshold=1.0e-6, + sync_type="alltoallv", + ) + cov_and_hits.apply(data) + + binner = ops.BinMap( + pixel_dist="pixel_dist", + covariance=cov_and_hits.covariance, + det_flags=None, + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + sync_type="alltoallv", + ) + + binner.det_data = key1 + binner.binned = "binned1" + binner.apply(data) + binner.det_data = key2 + binner.binned = "binned2" + binner.apply(data) + + # Study the map on the root process + + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + write_healpix_fits(data["binned1"], toast_bin_path, nest=pixels.nest) + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key2}.fits") + write_healpix_fits(data["binned2"], toast_bin_path, nest=pixels.nest) + + toast_hits_path = os.path.join(self.outdir, "toast_hits.fits") + write_healpix_fits(data[cov_and_hits.hits], toast_hits_path, nest=pixels.nest) + + fail = False + if self.rank == 0: + mapfile = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + mdata = hp.read_map(mapfile, field=range(3)) + mapfile = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + mdataw = hp.read_map(mapfile, field=range(3)) + + cl_out = hp.anafast(mdata, lmax=self.lmax) + cl_outw = hp.anafast(mdataw, lmax=self.lmax) + + np.testing.assert_almost_equal(cl_out, cl_outw, decimal=2) return - """ + def test_sim_TEB_conviqt(self): + if not ops.conviqt.available(): + print("libconviqt not available, skipping tests") + return + + # Create a fake scan strategy that hits every pixel once. + # data = create_healpix_ring_satellite(self.comm, nside=self.nside) + data = create_satellite_data( + self.comm, obs_time=120 * u.min, pixel_per_process=2 + ) + self.make_beam_file_dict(data) + + + # Generate timestreams + detpointing = ops.PointingDetectorSimple() + + + key1 = "conviqt0" + sim_conviqt = ops.SimConviqt( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False, + det_data=key1, + pol=True, + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + sim_conviqt.apply(data) + + key2 = "tebconviqt" + + sim_wconviqt = ops.SimTEBConviqt( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False , + det_data=key2, + pol=False , + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + + sim_wconviqt.apply(data) + # Bin a map to study + + pixels = ops.PixelsHealpix( + nside=self.nside, + nest=False, + detector_pointing=detpointing, + ) + pixels.apply(data) + weights = ops.StokesWeights( + mode="IQU", + hwp_angle=None, + detector_pointing=detpointing, + ) + weights.apply(data) + + default_model = ops.DefaultNoiseModel() + default_model.apply(data) + + cov_and_hits = ops.CovarianceAndHits( + pixel_dist="pixel_dist", + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + rcond_threshold=1.0e-6, + sync_type="alltoallv", + ) + cov_and_hits.apply(data) + + binner = ops.BinMap( + pixel_dist="pixel_dist", + covariance=cov_and_hits.covariance, + det_flags=None, + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + sync_type="alltoallv", + ) + + binner.det_data = key1 + binner.binned = "binned1" + binner.apply(data) + binner.det_data = key2 + binner.binned = "binned2" + binner.apply(data) + + # Study the map on the root process + + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + write_healpix_fits(data["binned1"], toast_bin_path, nest=pixels.nest) + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key2}.fits") + write_healpix_fits(data["binned2"], toast_bin_path, nest=pixels.nest) + + toast_hits_path = os.path.join(self.outdir, "toast_hits.fits") + write_healpix_fits(data[cov_and_hits.hits], toast_hits_path, nest=pixels.nest) + + fail = False + if self.rank == 0: + mapfile = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + mdata = hp.read_map(mapfile, field=range(3)) + mapfile = os.path.join(self.outdir, f"toast_bin.{key2}.fits") + mdataw = hp.read_map(mapfile, field=range(3)) + + cl_out = hp.anafast(mdata, lmax=self.lmax) + cl_outw = hp.anafast(mdataw, lmax=self.lmax) + + np.testing.assert_almost_equal(cl_out, cl_outw, decimal=2) + + return + + """ def test_sim_hwp(self): if not ops.conviqt.available(): print("libconviqt not available, skipping tests") return - # Create a fake scan strategy that hits every pixel once. data = create_satellite_data(self.comm) - # make a simple pointing matrix detpointing = ops.PointingDetectorSimple() - # Generate timestreams sim_conviqt = ops.SimWeightedConviqt( comm=self.comm, @@ -305,7 +587,5 @@ def test_sim_hwp(self): hwp_angle="hwp_angle", ) sim_conviqt.exec(data) - return - """ From 84287bf7e7fdd879f8a3bbe05b738de68e2a2fab Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 26 Apr 2022 14:17:06 +0200 Subject: [PATCH 2/3] port TEB changes to totalconvolve --- src/toast/ops/__init__.py | 6 +- src/toast/ops/totalconvolve.py | 239 +++++++ .../tests/ops_sim_tod_totalconvolve_new.py | 591 ++++++++++++++++++ 3 files changed, 835 insertions(+), 1 deletion(-) create mode 100644 src/toast/tests/ops_sim_tod_totalconvolve_new.py diff --git a/src/toast/ops/__init__.py b/src/toast/ops/__init__.py index d56c7d8f5..19af97324 100644 --- a/src/toast/ops/__init__.py +++ b/src/toast/ops/__init__.py @@ -68,7 +68,11 @@ from .sim_cosmic_rays import InjectCosmicRays -from .totalconvolve import SimTotalconvolve +from .totalconvolve import ( + SimTotalconvolve, + SimWeightedTotalconvolve, + SimTEBTotalconvolve, +) from .polyfilter import PolyFilter, PolyFilter2D, CommonModeFilter diff --git a/src/toast/ops/totalconvolve.py b/src/toast/ops/totalconvolve.py index 126296ff0..4cc9dfaec 100644 --- a/src/toast/ops/totalconvolve.py +++ b/src/toast/ops/totalconvolve.py @@ -277,6 +277,13 @@ def _exec(self, data, detectors=None, **kwargs): beam = self.get_beam(beam_file, lmax, mmax, det, verbose) theta, phi, psi, psi_pol = self.get_pointing(data, det, verbose) + np.savez( + f"{self.det_data}_{det}.npz", + psi=psi, + psi_pol=psi_pol, + theta=theta, + phi=phi, + ) pnt = self.get_buffer(theta, phi, psi, det, verbose) del theta, phi, psi if self.hwp_angle is None: @@ -820,3 +827,235 @@ def _provides(self): prov = self.detector_pointing.provides() prov["detdata"].append(self.det_data) return prov + + +class SimWeightedTotalconvolve(SimTotalconvolve): + """Operator which uses ducc0 to generate beam-convolved timestreams. + This operator should be used in presence of a spinning HWP which makes + the beam time-dependent, constantly mapping the co- and cross polar + responses on to each other. In OpSimTotalconvolve we assume the beam + to be static. + """ + + @function_timer + def _exec(self, data, detectors=None, **kwargs): + if not self.available: + raise RuntimeError("ducc0.totalconvolve is not available") + + if self.comm is None: + raise RuntimeError("ducc0.totalconvolve requires MPI") + + if self.detector_pointing is None: + raise RuntimeError("detector_pointing cannot be None.") + + log = Logger.get() + + timer = Timer() + timer.start() + + # Expand detector pointing + self.detector_pointing.apply(data, detectors=detectors) + + all_detectors = self._get_all_detectors(data, detectors) + + env = Environment.get() + nthreads = env.max_threads() + + for det in all_detectors: + verbose = self.comm.rank == 0 and self.verbosity > 0 + + # Expand detector pointing + self.detector_pointing.apply(data, detectors=[det]) + + if det in self.sky_file_dict: + sky_file = self.sky_file_dict[det] + else: + sky_file = self.sky_file.format(detector=det, mc=self.mc) + + if det in self.beam_file_dict: + beam_file = self.beam_file_dict[det] + else: + beam_file = self.beam_file.format(detector=det, mc=self.mc) + + lmax, mmax = self.get_lmmax( + sky_file, beam_file.replace(".fits", "_I000.fits") + ) + sky = self.get_sky(sky_file, lmax, det, verbose) + beamI00, beam0I0, beam00I = self.get_beam_weighted( + beam_file, lmax, mmax, det, verbose + ) + + theta, phi, psi, psi_pol = self.get_pointing(data, det, verbose) + + # I-beam convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data = self.convolve( + sky, beamI00, lmax, mmax, pnt, psi_pol, det, nthreads, verbose + ) + del pnt + + # Q-beam convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.cos(2 * psi_pol) * self.convolve( + sky, beam0I0, lmax, mmax, pnt, psi_pol, det, nthreads, verbose + ) + del pnt + + # U-beam convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.sin(2 * psi_pol) * self.convolve( + sky, beam00I, lmax, mmax, pnt, psi_pol, det, nthreads, verbose + ) + del theta, phi, psi + + self.calibrate_signal( + data, + det, + beamI00, + convolved_data, + verbose, + ) + self.save(data, det, convolved_data, verbose) + + del pnt, beamI00, beam0I0, beam00I, sky + + if verbose: + timer.report_clear(f"conviqt process detector {det}") + + return + + def get_beam_weighted(self, beamfile, lmax, mmax, det, verbose): + timer = Timer() + timer.start() + beam_file_i00 = beamfile.replace(".fits", "_I000.fits") + beam_file_0i0 = beamfile.replace(".fits", "_0I00.fits") + beam_file_00i = beamfile.replace(".fits", "_00I0.fits") + beami00 = self.get_beam(beam_file_i00, lmax, mmax, det, verbose) + beam0i0 = self.get_beam(beam_file_0i0, lmax, mmax, det, verbose) + beam00i = self.get_beam(beam_file_00i, lmax, mmax, det, verbose) + + if verbose: + timer.report_clear(f"initialize beam for detector {det}") + return beami00, beam0i0, beam00i + + +class SimTEBTotalconvolve(SimTotalconvolve): + """ + Operator that uses ducc0.totalconvolve to generate beam-convolved timestreams. + This operator should be used in presence of a spinning HWP which makes the beam time-dependent, + constantly mapping the co- and cross-polar responses on to each other. + In the parent class OpSimConviqt we assume the beam to be static. + + + The convolution is performed by coupling each IQU component of the signal propertly as: + :math:`skyT_lm * beamT_lm, skyE_lm * Re{P}, skyB_lm * Im{P}`. + + For extra details please refer to [this note ](https://giuspugl.github.io/reports/Notes_TEB_convolution.html) + """ + + @function_timer + def _exec(self, data, detectors=None, **kwargs): + if not self.available: + raise RuntimeError("ducc0.totalconvolve is not available") + + if self.comm is None: + raise RuntimeError("ducc0.totalconvolve requires MPI") + + if self.detector_pointing is None: + raise RuntimeError("detector_pointing cannot be None.") + + log = Logger.get() + + timer = Timer() + timer.start() + + # Expand detector pointing + self.detector_pointing.apply(data, detectors=detectors) + + all_detectors = self._get_all_detectors(data, detectors) + + env = Environment.get() + nthreads = env.max_threads() + + for det in all_detectors: + verbose = self.comm.rank == 0 and self.verbosity > 0 + + # Expand detector pointing + self.detector_pointing.apply(data, detectors=[det]) + + if det in self.sky_file_dict: + sky_file = self.sky_file_dict[det] + else: + sky_file = self.sky_file.format(detector=det, mc=self.mc) + + if det in self.beam_file_dict: + beam_file = self.beam_file_dict[det] + else: + beam_file = self.beam_file.format(detector=det, mc=self.mc) + + lmax, mmax = self.get_lmmax( + sky_file.replace(".fits", "_T.fits"), + beam_file.replace(".fits", "_T.fits"), + ) + + skyT = self.get_sky( + sky_file.replace(".fits", "_T.fits"), lmax, det, verbose + ) + skyE = self.get_sky( + sky_file.replace(".fits", "_E.fits"), lmax, det, verbose + ) + skyB = self.get_sky( + sky_file.replace(".fits", "_B.fits"), lmax, det, verbose + ) + + beam_T, beam_P = self.get_beam_TEB(beam_file, lmax, mmax, det, verbose) + + theta, phi, psi, psi_pol = self.get_pointing(data, det, verbose) + # T-convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + + convolved_data = self.convolve( + skyT, beam_T, lmax, mmax, pnt, psi_pol, det, nthreads, verbose + ) + + del (pnt,) + # E-convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.cos(2 * psi_pol) * self.convolve( + skyE, beam_P, lmax, mmax, pnt, psi_pol, det, nthreads, verbose + ) + del (pnt,) + # B-convolution + pnt = self.get_buffer(theta, phi, psi, det, verbose) + convolved_data += np.sin(2 * psi_pol) * self.convolve( + skyB, beam_P, lmax, mmax, pnt, psi_pol, det, nthreads, verbose + ) + del theta, phi, psi + + self.calibrate_signal( + data, + det, + beam_T, + convolved_data, + verbose, + ) + self.save(data, det, convolved_data, verbose) + + del pnt, beam_T, beam_P, skyT, skyE, skyB + + if verbose: + timer.report_clear(f"conviqt process detector {det}") + + return + + def get_beam_TEB(self, beamfile, lmax, mmax, det, verbose): + timer = Timer() + timer.start() + beam_file_T = beamfile.replace(".fits", "_T.fits") + beam_file_P = beamfile.replace(".fits", "_P.fits") + beamT = self.get_beam(beam_file_T, lmax, mmax, det, verbose) + beamP = self.get_beam(beam_file_P, lmax, mmax, det, verbose) + + if verbose: + timer.report_clear(f"initialize beam for detector {det}") + return beamT, beamP diff --git a/src/toast/tests/ops_sim_tod_totalconvolve_new.py b/src/toast/tests/ops_sim_tod_totalconvolve_new.py new file mode 100644 index 000000000..658e72d85 --- /dev/null +++ b/src/toast/tests/ops_sim_tod_totalconvolve_new.py @@ -0,0 +1,591 @@ +# Copyright (c) 2015-2020 by the parties listed in the AUTHORS file. +# All rights reserved. Use of this source code is governed by +# a BSD-style license that can be found in the LICENSE file. + +import os + +import numpy as np + +from astropy import units as u + +import healpy as hp + +from .mpi import MPITestCase + +from ..vis import set_matplotlib_backend + +from .. import qarray as qa + +from .. import ops as ops + +from ..observation import default_values as defaults + +from ..pixels_io import write_healpix_fits + +from ._helpers import ( + create_outdir, + create_healpix_ring_satellite, + create_satellite_data, + create_fake_sky_alm, + create_fake_beam_alm, +) + + +class SimTotalconvolveTest(MPITestCase): + def setUp(self): + + np.random.seed(777) + fixture_name = os.path.splitext(os.path.basename(__file__))[0] + self.outdir = create_outdir(self.comm, fixture_name) + + self.nside = 32 + self.lmax = 128 + self.fwhm_sky = 10 * u.degree + self.fwhm_beam = 15 * u.degree + self.mmax = self.lmax + self.fname_sky = os.path.join(self.outdir, "sky_alm.fits") + self.fname_beam = os.path.join(self.outdir, "beam_alm.fits") + + self.rank = 0 + if self.comm is not None: + self.rank = self.comm.rank + + if self.rank == 0: + # Synthetic sky and beam (a_lm expansions) + self.slm = create_fake_sky_alm(self.lmax, self.fwhm_sky) + + hp.write_alm(self.fname_sky, self.slm, lmax=self.lmax, overwrite=True) + hp.write_alm(self.fname_sky.replace(".fits","_T.fits"), self.slm[0] , lmax=self.lmax, overwrite=True) + hp.write_alm(self.fname_sky.replace(".fits","_E.fits"), self.slm[1] , lmax=self.lmax, overwrite=True) + hp.write_alm(self.fname_sky.replace(".fits","_B.fits"), self.slm[2] , lmax=self.lmax, overwrite=True) + + self.blm = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + normalize_beam=True, + ) + + self.blm_bottom = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + normalize_beam=True, + detB_beam=True, + ) + + hp.write_alm( + self.fname_beam, + self.blm, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom.fits"), + self.blm_bottom, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + blm_I, blm_Q, blm_U = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + separate_IQU=True, + normalize_beam=True, + ) + blm_Ibot, blm_Qbot, blm_Ubot = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + separate_IQU=True, + detB_beam=True, + normalize_beam=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_I000.fits"), + blm_I, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_0I00.fits"), + blm_Q, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_00I0.fits"), + blm_U, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_I000.fits"), + blm_Ibot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_0I00.fits"), + blm_Qbot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_00I0.fits"), + blm_Ubot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + + # we explicitly store 3 separate beams for the T, E and B sky alm. + blm_T, blm_P = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + separate_TP=True, + detB_beam=False, + normalize_beam=True, + ) + blm_Tbot, blm_Pbot = create_fake_beam_alm( + self.lmax, + self.mmax, + fwhm_x=self.fwhm_beam, + fwhm_y=self.fwhm_beam, + separate_TP=True, + detB_beam=True, + normalize_beam=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_T.fits"), + blm_T, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_T.fits"), + blm_Tbot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_P.fits"), + blm_P, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + hp.write_alm( + self.fname_beam.replace(".fits", "_bottom_P.fits"), + blm_Pbot, + lmax=self.lmax, + mmax_in=self.mmax, + overwrite=True, + ) + # in + + if self.comm is not None: + self.comm.barrier() + + return + + def make_beam_file_dict(self, data): + """ + We make sure that data observed in each A/B detector within a + detector pair will have the right beam associated to it. In particular, + Q and U beams of B detector must have a flipped sign wrt the A one. + """ + + fname2 = self.fname_beam.replace(".fits", "_bottom.fits") + + self.beam_file_dict = {} + for det in data.obs[0].local_detectors: + if det[-1] == "A": + self.beam_file_dict[det] = self.fname_beam + else: + self.beam_file_dict[det] = fname2 + + return + + def test_sim_totalconvolve(self): + if not ops.totalconvolve.available(): + print("libtotalconvolve not available, skipping tests") + return + + # Create a fake scan strategy that hits every pixel once. + # data = create_healpix_ring_satellite(self.comm, nside=self.nside) + data = create_satellite_data( + self.comm, obs_time=120 * u.min, pixel_per_process=2 + ) + self.make_beam_file_dict(data) + + # Generate timestreams + + detpointing = ops.PointingDetectorSimple() + + key = defaults.det_data + sim_totalconvolve = ops.SimTotalconvolve( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False, + det_data=key, + pol=True, + normalize_beam=False, # beams are already produced normalized + fwhm=self.fwhm_sky, + ) + + sim_totalconvolve.apply(data) + + # Bin a map to study + + pixels = ops.PixelsHealpix( + nside=self.nside, + nest=False, + detector_pointing=detpointing, + ) + pixels.apply(data) + weights = ops.StokesWeights( + mode="IQU", + hwp_angle=None, + detector_pointing=detpointing, + ) + weights.apply(data) + + default_model = ops.DefaultNoiseModel() + default_model.apply(data) + + cov_and_hits = ops.CovarianceAndHits( + pixel_dist="pixel_dist", + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + rcond_threshold=1.0e-6, + sync_type="alltoallv", + ) + cov_and_hits.apply(data) + + binner = ops.BinMap( + pixel_dist="pixel_dist", + covariance=cov_and_hits.covariance, + det_data=key, + det_flags=None, + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + sync_type="alltoallv", + ) + binner.apply(data) + + # Study the map on the root process + + toast_bin_path = os.path.join(self.outdir, "toast_bin.fits") + write_healpix_fits(data[binner.binned], toast_bin_path, nest=pixels.nest) + + toast_hits_path = os.path.join(self.outdir, "toast_hits.fits") + write_healpix_fits(data[cov_and_hits.hits], toast_hits_path, nest=pixels.nest) + + fail = False + + if self.rank == 0: + hitsfile = os.path.join(self.outdir, "toast_hits.fits") + + hdata = hp.read_map(hitsfile) + + footprint = np.ma.masked_not_equal(hdata, 0.0).mask + + mapfile = os.path.join(self.outdir, "toast_bin.fits") + mdata = hp.read_map(mapfile, field=range(3)) + + deconv = 1 / hp.gauss_beam( + self.fwhm_sky.to_value(u.radian), + lmax=self.lmax, + pol=False, + ) + + smoothed = hp.alm2map( + [hp.almxfl(self.slm[ii], deconv) for ii in range(3)], + self.nside, + lmax=self.lmax, + fwhm=self.fwhm_beam.to_value(u.radian), + verbose=False, + pixwin=False, + ) + smoothed[:, ~footprint] = 0 + cl_out = hp.anafast(mdata, lmax=self.lmax) + cl_smoothed = hp.anafast(smoothed, lmax=self.lmax) + + np.testing.assert_almost_equal(cl_smoothed[0], cl_out[0], decimal=2) + + return + + def test_sim_weighted_totalconvolve(self): + if not ops.totalconvolve.available(): + print("libtotalconvolve not available, skipping tests") + return + + # Create a fake scan strategy that hits every pixel once. + # data = create_healpix_ring_satellite(self.comm, nside=self.nside) + data = create_satellite_data( + self.comm, obs_time=120 * u.min, pixel_per_process=2 + ) + self.make_beam_file_dict(data) + + # Generate timestreams + + detpointing = ops.PointingDetectorSimple() + + key1 = "totalconvolve" + sim_totalconvolve = ops.SimTotalconvolve( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False, + det_data=key1, + pol=True, + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + sim_totalconvolve.apply(data) + + key2 = "wtotalconvolve" + + sim_wtotalconvolve = ops.SimWeightedTotalconvolve( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False, + det_data=key2, + pol=True, + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + sim_wtotalconvolve.apply(data) + # Bin a map to study + + pixels = ops.PixelsHealpix( + nside=self.nside, + nest=False, + detector_pointing=detpointing, + ) + pixels.apply(data) + weights = ops.StokesWeights( + mode="IQU", + hwp_angle=None, + detector_pointing=detpointing, + ) + weights.apply(data) + + default_model = ops.DefaultNoiseModel() + default_model.apply(data) + + cov_and_hits = ops.CovarianceAndHits( + pixel_dist="pixel_dist", + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + rcond_threshold=1.0e-6, + sync_type="alltoallv", + ) + cov_and_hits.apply(data) + + binner = ops.BinMap( + pixel_dist="pixel_dist", + covariance=cov_and_hits.covariance, + det_flags=None, + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + sync_type="alltoallv", + ) + + binner.det_data = key1 + binner.binned = "binned1" + binner.apply(data) + binner.det_data = key2 + binner.binned = "binned2" + binner.apply(data) + + # Study the map on the root process + + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + write_healpix_fits(data["binned1"], toast_bin_path, nest=pixels.nest) + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key2}.fits") + write_healpix_fits(data["binned2"], toast_bin_path, nest=pixels.nest) + + toast_hits_path = os.path.join(self.outdir, "toast_hits.fits") + write_healpix_fits(data[cov_and_hits.hits], toast_hits_path, nest=pixels.nest) + + fail = False + if self.rank == 0: + mapfile = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + mdata = hp.read_map(mapfile, field=range(3)) + mapfile = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + mdataw = hp.read_map(mapfile, field=range(3)) + + cl_out = hp.anafast(mdata, lmax=self.lmax) + cl_outw = hp.anafast(mdataw, lmax=self.lmax) + + np.testing.assert_almost_equal(cl_out, cl_outw, decimal=2) + + return + + def test_sim_TEB_totalconvolve(self): + if not ops.totalconvolve.available(): + print("libtotalconvolve not available, skipping tests") + return + + # Create a fake scan strategy that hits every pixel once. + # data = create_healpix_ring_satellite(self.comm, nside=self.nside) + data = create_satellite_data( + self.comm, obs_time=120 * u.min, pixel_per_process=2 + ) + self.make_beam_file_dict(data) + + + # Generate timestreams + + detpointing = ops.PointingDetectorSimple() + + + key1 = "totalconvolve0" + sim_totalconvolve = ops.SimTotalconvolve( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False, + det_data=key1, + pol=True, + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + sim_totalconvolve.apply(data) + + key2 = "tebtotalconvolve" + + sim_wtotalconvolve = ops.SimTEBTotalconvolve( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file_dict=self.beam_file_dict, + dxx=False , + det_data=key2, + pol=False , + normalize_beam=False, + fwhm=self.fwhm_sky, + ) + + + sim_wtotalconvolve.apply(data) + # Bin a map to study + + pixels = ops.PixelsHealpix( + nside=self.nside, + nest=False, + detector_pointing=detpointing, + ) + pixels.apply(data) + weights = ops.StokesWeights( + mode="IQU", + hwp_angle=None, + detector_pointing=detpointing, + ) + weights.apply(data) + + default_model = ops.DefaultNoiseModel() + default_model.apply(data) + + cov_and_hits = ops.CovarianceAndHits( + pixel_dist="pixel_dist", + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + rcond_threshold=1.0e-6, + sync_type="alltoallv", + ) + cov_and_hits.apply(data) + + binner = ops.BinMap( + pixel_dist="pixel_dist", + covariance=cov_and_hits.covariance, + det_flags=None, + pixel_pointing=pixels, + stokes_weights=weights, + noise_model=default_model.noise_model, + sync_type="alltoallv", + ) + + binner.det_data = key1 + binner.binned = "binned1" + binner.apply(data) + binner.det_data = key2 + binner.binned = "binned2" + binner.apply(data) + + # Study the map on the root process + + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + write_healpix_fits(data["binned1"], toast_bin_path, nest=pixels.nest) + toast_bin_path = os.path.join(self.outdir, f"toast_bin.{key2}.fits") + write_healpix_fits(data["binned2"], toast_bin_path, nest=pixels.nest) + + toast_hits_path = os.path.join(self.outdir, "toast_hits.fits") + write_healpix_fits(data[cov_and_hits.hits], toast_hits_path, nest=pixels.nest) + + fail = False + if self.rank == 0: + mapfile = os.path.join(self.outdir, f"toast_bin.{key1}.fits") + mdata = hp.read_map(mapfile, field=range(3)) + mapfile = os.path.join(self.outdir, f"toast_bin.{key2}.fits") + mdataw = hp.read_map(mapfile, field=range(3)) + + cl_out = hp.anafast(mdata, lmax=self.lmax) + cl_outw = hp.anafast(mdataw, lmax=self.lmax) + + np.testing.assert_almost_equal(cl_out, cl_outw, decimal=2) + + return + + """ + def test_sim_hwp(self): + if not ops.totalconvolve.available(): + print("libtotalconvolve not available, skipping tests") + return + # Create a fake scan strategy that hits every pixel once. + data = create_satellite_data(self.comm) + # make a simple pointing matrix + detpointing = ops.PointingDetectorSimple() + # Generate timestreams + sim_totalconvolve = ops.SimWeightedTotalconvolve( + comm=self.comm, + detector_pointing=detpointing, + sky_file=self.fname_sky, + beam_file=self.fname_beam, + dxx=False, + hwp_angle="hwp_angle", + ) + sim_totalconvolve.exec(data) + return + """ From 49291f77cc6d5f7b0d859f8dca4d0b05d5ccf25b Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Thu, 28 Apr 2022 08:14:11 +0200 Subject: [PATCH 3/3] comment save statement --- src/toast/ops/totalconvolve.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/toast/ops/totalconvolve.py b/src/toast/ops/totalconvolve.py index 4cc9dfaec..92707a9e5 100644 --- a/src/toast/ops/totalconvolve.py +++ b/src/toast/ops/totalconvolve.py @@ -277,13 +277,13 @@ def _exec(self, data, detectors=None, **kwargs): beam = self.get_beam(beam_file, lmax, mmax, det, verbose) theta, phi, psi, psi_pol = self.get_pointing(data, det, verbose) - np.savez( - f"{self.det_data}_{det}.npz", - psi=psi, - psi_pol=psi_pol, - theta=theta, - phi=phi, - ) +# np.savez( +# f"{self.det_data}_{det}.npz", +# psi=psi, +# psi_pol=psi_pol, +# theta=theta, +# phi=phi, +# ) pnt = self.get_buffer(theta, phi, psi, det, verbose) del theta, phi, psi if self.hwp_angle is None: