diff --git a/.gitignore b/.gitignore index e35fea85..5eeaa6cd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,29 @@ +### PyCharm ### +# Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio, WebStorm and Rider +# Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839 + +# User-specific stuff +.idea/**/workspace.xml +.idea/**/tasks.xml +.idea/**/usage.statistics.xml +.idea/**/dictionaries +.idea/**/shelf + +# AWS User-specific +.idea/**/aws.xml + +# Generated files +.idea/**/contentModel.xml + +# Sensitive or high-churn files +.idea/**/dataSources/ +.idea/**/dataSources.ids +.idea/**/dataSources.local.xml +.idea/**/sqlDataSources.xml +.idea/**/dynamic.xml +.idea/**/uiDesigner.xml +.idea/**/dbnavigator.xml + *~ # Byte-compiled / optimized / DLL files diff --git a/.idea/inspectionProfiles/Project_Default.xml b/.idea/inspectionProfiles/Project_Default.xml new file mode 100644 index 00000000..ad2dd076 --- /dev/null +++ b/.idea/inspectionProfiles/Project_Default.xml @@ -0,0 +1,17 @@ + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 00000000..105ce2da --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/litebird_sim.iml b/.idea/litebird_sim.iml new file mode 100644 index 00000000..07f841d1 --- /dev/null +++ b/.idea/litebird_sim.iml @@ -0,0 +1,15 @@ + + + + + + + + + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 00000000..c26b3564 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,7 @@ + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 00000000..f6c88267 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 00000000..94a25f7f --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index 38f27d56..32db5dc9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,8 @@ - Each `Simulation` object creates random number generators (field `Simulation.random`), in a way that is safe even for MPI applications +- Support the production of maps in Galactic coordinates through the TOAST2 map-maker + - Make `make_bin_map` compute pixel indices instead of requiring them as input [#176](https://github.com/litebird/litebird_sim/pull/176) - Use a more robust algorithm to compute pointings [#175](https://github.com/litebird/litebird_sim/pull/175) diff --git a/litebird_sim/destriper/__init__.py b/litebird_sim/destriper/__init__.py index 1b02962b..263dcffe 100644 --- a/litebird_sim/destriper/__init__.py +++ b/litebird_sim/destriper/__init__.py @@ -15,6 +15,8 @@ from toast.tod.interval import Interval import toast.mpi +from litebird_sim.coordinates import CoordinateSystem, rotate_coordinates_e2g + toast.mpi.use_mpi = lbs.MPI_ENABLED @@ -26,6 +28,10 @@ class DestriperParameters: - ``nside``: the NSIDE parameter used to create the maps + - ``coordinate_system``: an instance of the :class:`.CoordinateSystem` enum. + It specifies if the map must be created in ecliptic (default) or + galactic coordinates. + - ``nnz``: number of components per pixel. The default is 3 (I/Q/U). - ``baseline_length``: number of consecutive samples in a 1/f noise @@ -64,6 +70,7 @@ class DestriperParameters: """ nside: int = 512 + coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic nnz: int = 3 baseline_length: int = 100 iter_max: int = 100 @@ -104,12 +111,20 @@ class DestriperResult: npp: Any = None invnpp: Any = None rcond: Any = None + coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic class _Toast2FakeCache: "This class simulates a TOAST2 cache" - def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside): + def __init__( + self, + spin2ecliptic_quats, + obs, + bore2spin_quat, + nside, + coordinates: CoordinateSystem, + ): self.obs = obs self.keydict = {"timestamps": obs.get_times()} @@ -141,13 +156,27 @@ def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside): logging.warning( "converting TODs for %s from %s to float64", obs.name[i], - str(pointings[i].dtype), + str(obs.tod[i].dtype), ) self.keydict[f"signal_{det}"] = np.array(obs.tod[i], dtype=np.float64) - self.keydict[f"pixels_{det}"] = healpix_base.ang2pix(curpnt[:, 0:2]) + theta_phi = curpnt[:, 0:2] + polangle = curpnt[:, 2] + + if coordinates == CoordinateSystem.Galactic: + theta_phi, polangle = rotate_coordinates_e2g( + pointings_ecl=theta_phi, pol_angle_ecl=polangle + ) + elif coordinates == CoordinateSystem.Ecliptic: + pass # Do nothing, "theta_phi" and "polangle" are ok + else: + assert ValueError( + "unable to handle coordinate system {coordinates} in `destripe`" + ) + + self.keydict[f"pixels_{det}"] = healpix_base.ang2pix(theta_phi) self.keydict[f"weights_{det}"] = np.stack( - (np.ones(nsamples), np.cos(2 * curpnt[:, 2]), np.sin(2 * curpnt[:, 2])) + (np.ones(nsamples), np.cos(2 * polangle), np.sin(2 * polangle)) ).transpose() def keys(self): @@ -172,10 +201,19 @@ def destroy(self, name): class _Toast2FakeTod: "This class simulates a TOAST2 TOD" - def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside): + def __init__( + self, + spin2ecliptic_quats, + obs, + bore2spin_quat, + nside, + coordinates: CoordinateSystem, + ): self.obs = obs self.local_samples = (0, obs.tod[0].size) - self.cache = _Toast2FakeCache(spin2ecliptic_quats, obs, bore2spin_quat, nside) + self.cache = _Toast2FakeCache( + spin2ecliptic_quats, obs, bore2spin_quat, nside, coordinates + ) def local_intervals(self, _): return [ @@ -215,9 +253,20 @@ def local_times(self): class _Toast2FakeData: "This class simulates a TOAST2 Data class" - def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside): + def __init__( + self, + spin2ecliptic_quats, + obs, + bore2spin_quat, + nside, + coordinates: CoordinateSystem, + ): self.obs = [ - {"tod": _Toast2FakeTod(spin2ecliptic_quats, x, bore2spin_quat, nside)} + { + "tod": _Toast2FakeTod( + spin2ecliptic_quats, x, bore2spin_quat, nside, coordinates + ) + } for x in obs ] self.bore2spin_quat = bore2spin_quat @@ -285,6 +334,7 @@ def destripe_observations( obs=observations, bore2spin_quat=bore2spin_quat, nside=params.nside, + coordinates=params.coordinate_system, ) mapmaker = OpMapMaker( nside=params.nside, @@ -306,6 +356,8 @@ def destripe_observations( result = DestriperResult() + result.coordinate_system = params.coordinate_system + if params.return_hit_map: result.hit_map = healpy.read_map( base_path / (params.output_file_prefix + "hits.fits"), diff --git a/litebird_sim/scan_map.py b/litebird_sim/scan_map.py index 9ec206c5..70771286 100644 --- a/litebird_sim/scan_map.py +++ b/litebird_sim/scan_map.py @@ -140,7 +140,7 @@ def scan_map_in_observations( elif all(item in maps.keys() for item in cur_obs.channel): input_names = cur_obs.channel else: - print( + raise ValueError( "The dictionary maps does not contain all the relevant" + "keys, please check the list of detectors and channels" ) diff --git a/test/destriper_reference/destriper_binned_map.fits.gz b/test/destriper_reference/destriper_binned_map_ecliptic.fits.gz similarity index 100% rename from test/destriper_reference/destriper_binned_map.fits.gz rename to test/destriper_reference/destriper_binned_map_ecliptic.fits.gz diff --git a/test/destriper_reference/destriper_binned_map_galactic.fits.gz b/test/destriper_reference/destriper_binned_map_galactic.fits.gz new file mode 100644 index 00000000..3154bbd1 Binary files /dev/null and b/test/destriper_reference/destriper_binned_map_galactic.fits.gz differ diff --git a/test/destriper_reference/destriper_destriped_map.fits.gz b/test/destriper_reference/destriper_destriped_map_ecliptic.fits.gz similarity index 100% rename from test/destriper_reference/destriper_destriped_map.fits.gz rename to test/destriper_reference/destriper_destriped_map_ecliptic.fits.gz diff --git a/test/destriper_reference/destriper_destriped_map_galactic.fits.gz b/test/destriper_reference/destriper_destriped_map_galactic.fits.gz new file mode 100644 index 00000000..48cb4480 Binary files /dev/null and b/test/destriper_reference/destriper_destriped_map_galactic.fits.gz differ diff --git a/test/destriper_reference/destriper_hit_map.fits.gz b/test/destriper_reference/destriper_hit_map_ecliptic.fits.gz similarity index 100% rename from test/destriper_reference/destriper_hit_map.fits.gz rename to test/destriper_reference/destriper_hit_map_ecliptic.fits.gz diff --git a/test/destriper_reference/destriper_hit_map_galactic.fits.gz b/test/destriper_reference/destriper_hit_map_galactic.fits.gz new file mode 100644 index 00000000..091eb579 Binary files /dev/null and b/test/destriper_reference/destriper_hit_map_galactic.fits.gz differ diff --git a/test/destriper_reference/destriper_invnpp.fits.gz b/test/destriper_reference/destriper_invnpp_ecliptic.fits.gz similarity index 100% rename from test/destriper_reference/destriper_invnpp.fits.gz rename to test/destriper_reference/destriper_invnpp_ecliptic.fits.gz diff --git a/test/destriper_reference/destriper_invnpp_galactic.fits.gz b/test/destriper_reference/destriper_invnpp_galactic.fits.gz new file mode 100644 index 00000000..c73a07d5 Binary files /dev/null and b/test/destriper_reference/destriper_invnpp_galactic.fits.gz differ diff --git a/test/destriper_reference/destriper_npp.fits.gz b/test/destriper_reference/destriper_npp_ecliptic.fits.gz similarity index 100% rename from test/destriper_reference/destriper_npp.fits.gz rename to test/destriper_reference/destriper_npp_ecliptic.fits.gz diff --git a/test/destriper_reference/destriper_npp_galactic.fits.gz b/test/destriper_reference/destriper_npp_galactic.fits.gz new file mode 100644 index 00000000..d77aa5b6 Binary files /dev/null and b/test/destriper_reference/destriper_npp_galactic.fits.gz differ diff --git a/test/destriper_reference/destriper_rcond.fits.gz b/test/destriper_reference/destriper_rcond_ecliptic.fits.gz similarity index 100% rename from test/destriper_reference/destriper_rcond.fits.gz rename to test/destriper_reference/destriper_rcond_ecliptic.fits.gz diff --git a/test/destriper_reference/destriper_rcond_galactic.fits.gz b/test/destriper_reference/destriper_rcond_galactic.fits.gz new file mode 100644 index 00000000..8c09f0a2 Binary files /dev/null and b/test/destriper_reference/destriper_rcond_galactic.fits.gz differ diff --git a/test/test_destriper.py b/test/test_destriper.py index 866eca80..876aafdc 100644 --- a/test/test_destriper.py +++ b/test/test_destriper.py @@ -9,8 +9,18 @@ import healpy from numpy.random import MT19937, RandomState, SeedSequence +from litebird_sim import CoordinateSystem + + +COORDINATE_SYSTEM_STR = { + CoordinateSystem.Ecliptic: "ecliptic", + CoordinateSystem.Galactic: "galactic", +} + + +def run_destriper_tests(tmp_path, coordinates: CoordinateSystem): + coordinates_str = COORDINATE_SYSTEM_STR[coordinates] -def test_destriper(tmp_path): sim = lbs.Simulation( base_path=tmp_path / "destriper_output", start_time=0, duration_s=86400.0 ) @@ -49,6 +59,7 @@ def test_destriper(tmp_path): nnz=3, baseline_length=100, iter_max=10, + coordinate_system=coordinates, return_hit_map=True, return_binned_map=True, return_destriped_map=True, @@ -58,16 +69,19 @@ def test_destriper(tmp_path): ) results = lbs.destripe(sim, instr, params=params) + assert results.coordinate_system == coordinates ref_map_path = Path(__file__).parent / "destriper_reference" - hit_map_filename = ref_map_path / "destriper_hit_map.fits.gz" + hit_map_filename = ref_map_path / f"destriper_hit_map_{coordinates_str}.fits.gz" # healpy.write_map(hit_map_filename, results.hit_map, dtype="int32", overwrite=True) np.testing.assert_allclose( results.hit_map, healpy.read_map(hit_map_filename, field=None, dtype=np.int32) ) - binned_map_filename = ref_map_path / "destriper_binned_map.fits.gz" + binned_map_filename = ( + ref_map_path / f"destriper_binned_map_{coordinates_str}.fits.gz" + ) # healpy.write_map( # binned_map_filename, # results.binned_map, @@ -80,7 +94,9 @@ def test_destriper(tmp_path): assert results.binned_map.shape == ref_binned.shape np.testing.assert_allclose(results.binned_map, ref_binned, rtol=1e-2, atol=1e-3) - destriped_map_filename = ref_map_path / "destriper_destriped_map.fits.gz" + destriped_map_filename = ( + ref_map_path / f"destriper_destriped_map_{coordinates_str}.fits.gz" + ) # healpy.write_map( # destriped_map_filename, # results.destriped_map, @@ -95,7 +111,7 @@ def test_destriper(tmp_path): results.destriped_map, ref_destriped, rtol=1e-2, atol=1e-3 ) - npp_filename = ref_map_path / "destriper_npp.fits.gz" + npp_filename = ref_map_path / f"destriper_npp_{coordinates_str}.fits.gz" # healpy.write_map( # npp_filename, # results.npp, @@ -108,7 +124,7 @@ def test_destriper(tmp_path): assert results.npp.shape == ref_npp.shape np.testing.assert_allclose(results.npp, ref_npp, rtol=1e-2, atol=1e-3) - invnpp_filename = ref_map_path / "destriper_invnpp.fits.gz" + invnpp_filename = ref_map_path / f"destriper_invnpp_{coordinates_str}.fits.gz" # healpy.write_map( # invnpp_filename, # results.invnpp, @@ -121,7 +137,7 @@ def test_destriper(tmp_path): assert results.invnpp.shape == ref_invnpp.shape np.testing.assert_allclose(results.invnpp, ref_invnpp, rtol=1e-2, atol=1e-3) - rcond_filename = ref_map_path / "destriper_rcond.fits.gz" + rcond_filename = ref_map_path / f"destriper_rcond_{coordinates_str}.fits.gz" # healpy.write_map( # rcond_filename, # results.rcond, @@ -131,3 +147,112 @@ def test_destriper(tmp_path): assert np.allclose( results.rcond, healpy.read_map(rcond_filename, field=None, dtype=np.float32) ) + + +def test_destriper_ecliptic(tmp_path): + run_destriper_tests(tmp_path=tmp_path, coordinates=CoordinateSystem.Ecliptic) + + +def test_destriper_galactic(tmp_path): + run_destriper_tests(tmp_path=tmp_path, coordinates=CoordinateSystem.Galactic) + + +def test_destriper_coordinate_consistency(tmp_path): + # Here we check that MBS uses the same coordinate system as the destriper + # in «Galactic» mode: specifically, we create a noiseless TOD from a CMB + # map in Galactic coordinates and run the destriper asking to use Galactic + # coordinates again. Since the TOD was noiseless, the binned map should be + # the same as the input map, except for two features: + # + # 1. It does not cover the whole sky + # 2. A few pixels at the border of the observed region are not properly + # constrained and their value will be set to zero + # + # This test checks that «most» of the two maps agree, i.e., the 5% and 95% + # percentiles are smaller than some (small) threshold. + + sim = lbs.Simulation( + base_path="destriper_output", + start_time=0, + duration_s=10 * 86400.0, + ) + + sim.generate_spin2ecl_quaternions( + scanning_strategy=lbs.SpinningScanningStrategy( + spin_sun_angle_rad=np.deg2rad(30), # CORE-specific parameter + spin_rate_hz=0.5 / 60, # Ditto + # We use astropy to convert the period (4 days) in + # seconds + precession_rate_hz=1.0 / (4 * u.day).to("s").value, + ) + ) + instr = lbs.InstrumentInfo( + name="mock_instrument", + spin_boresight_angle_rad=np.deg2rad(65), + ) + + detectors = [ + lbs.DetectorInfo( + name="noiseless_detector", + sampling_rate_hz=10.0, + fwhm_arcmin=60.0, + bandcenter_ghz=40.0, + bandwidth_ghz=12.0, + ), + ] + + # We create two detectors, whose polarization angles are separated by π/2 + sim.create_observations( + detectors=detectors, + dtype_tod=np.float64, # Needed if you use the TOAST destriper + n_blocks_time=lbs.MPI_COMM_WORLD.size, + split_list_over_processes=False, + ) + + for obs in sim.observations: + lbs.get_pointings( + obs, + spin2ecliptic_quats=sim.spin2ecliptic_quats, + detector_quats=None, + bore2spin_quat=instr.bore2spin_quat, + ) + + params = lbs.MbsParameters( + make_cmb=True, + ) + mbs = lbs.Mbs( + simulation=sim, + parameters=params, + detector_list=detectors, + ) + (healpix_maps, file_paths) = mbs.run_all() + + lbs.scan_map_in_observations(obs=sim.observations, maps=healpix_maps) + + params = lbs.DestriperParameters( + nside=healpy.npix2nside(len(healpix_maps[detectors[0].name][0])), + coordinate_system=lbs.CoordinateSystem.Galactic, + return_hit_map=True, + return_binned_map=True, + return_destriped_map=True, + ) + + result = lbs.destripe(sim, instr, params) + + inp = healpix_maps[detectors[0].name] # Input CMB map in Galactic coordinates + out = result.binned_map # The binned map produced by the destriper + hit = result.hit_map # The hit map produced by the destriper + + # We do not consider unseen pixels nor pixels that have not been properly + # constrained by our mock detector + mask = (hit == 0) | (out == 0.0) + inp[mask] = np.NaN + out[mask] = np.NaN + + # Ideally this should be ≈0 + diff = inp - out + + # We check the closeness of `diff` to zero through the 5% and 95% percentiles + low_limit, high_limit = np.percentile(diff[np.isfinite(diff)], (0.05, 0.95)) + assert np.abs(low_limit) < 1e-9 + assert np.abs(high_limit) < 1e-9