From f4a097d89fb4fb0ba50deb377b23cc867b29f359 Mon Sep 17 00:00:00 2001 From: Nicolas Tessore Date: Sun, 12 Nov 2023 14:33:41 +0000 Subject: [PATCH 1/2] add "label" property to catalogs --- heracles/catalog/base.py | 21 +++++++++++++++++++++ tests/test_catalog.py | 9 +++++++++ 2 files changed, 30 insertions(+) diff --git a/heracles/catalog/base.py b/heracles/catalog/base.py index 4a79bf7..f505e5b 100644 --- a/heracles/catalog/base.py +++ b/heracles/catalog/base.py @@ -117,6 +117,11 @@ def __getitem__(self, where): """create a view with the given selection""" ... + @property + def label(self): + """return a human-friendly identifier for the source of data""" + ... + @property def base(self): """return the base catalogue of a view, or ``None`` if not a view""" @@ -186,6 +191,11 @@ def base(self): """base catalogue of this view""" return self._catalog + @property + def label(self): + """human-friendly label of the catalogue (not settable in view)""" + return self._catalog.label + @property def selection(self): """selection of this view""" @@ -251,6 +261,7 @@ def __init__(self): self._page_size = self.default_page_size self._filters = [] + self._label = None self._visibility = None def __copy__(self): @@ -259,6 +270,7 @@ def __copy__(self): other = self.__class__.__new__(self.__class__) other._page_size = self._page_size other._filters = self._filters.copy() + other._label = self._label other._visibility = self._visibility return other @@ -304,6 +316,15 @@ def base(self): """returns ``None`` since this is not a view of another catalogue""" return + @property + def label(self): + """optional human-friendly label for catalogue""" + return self._label + + @label.setter + def label(self, label): + self._label = label + @property def selection(self): """returns ``None`` since this is not a view of another catalogue""" diff --git a/tests/test_catalog.py b/tests/test_catalog.py index 8ab0ba9..55a9e76 100644 --- a/tests/test_catalog.py +++ b/tests/test_catalog.py @@ -180,6 +180,10 @@ def test_catalog_base_properties(catalog): catalog.filters = [] assert catalog.filters == [] + assert catalog.label is None + catalog.label = "label 123" + assert catalog.label == "label 123" + v = object() assert catalog.visibility is None catalog.visibility = v @@ -235,6 +239,7 @@ def _pages(self, selection): def test_catalog_view(catalog): from heracles.catalog import Catalog + catalog.label = "label 123" catalog.visibility = cvis = object() where = object() @@ -247,9 +252,13 @@ def test_catalog_view(catalog): assert catalog.base is None assert catalog.selection is None assert view.base is catalog + assert view.label == "label 123" assert view.selection is where assert view.visibility is catalog.visibility + with pytest.raises(AttributeError): + view.label = "different label" + view.visibility = vvis = object() assert view.visibility is not catalog.visibility From 0154ef56df4c0e98296f51a2b1655342be4c4a15 Mon Sep 17 00:00:00 2001 From: Nicolas Tessore Date: Sun, 12 Nov 2023 15:07:46 +0000 Subject: [PATCH 2/2] store input catalogue label in map metadata --- heracles/io.py | 6 ++ heracles/maps.py | 29 ++++++-- tests/test_io.py | 16 +++-- tests/test_maps.py | 160 ++++++++++++++++++++++++--------------------- 4 files changed, 127 insertions(+), 84 deletions(-) diff --git a/heracles/io.py b/heracles/io.py index a40c0f5..aad3fa2 100644 --- a/heracles/io.py +++ b/heracles/io.py @@ -36,14 +36,20 @@ _METADATA_COMMENTS = { + "catalog": "catalog of map", "spin": "spin weight of map", "kernel": "mapping kernel of map", + "nside": "NSIDE parameter of HEALPix map", "power": "area power of map", + "catalog_1": "catalog of first map", "spin_1": "spin weight of first map", "kernel_1": "mapping kernel of first map", + "nside_1": "NSIDE parameter of first HEALPix map", "power_1": "area power of first map", + "catalog_2": "catalog of second map", "spin_2": "spin weight of second map", "kernel_2": "mapping kernel of second map", + "nside_2": "NSIDE parameter of second HEALPix map", "power_2": "area power of second map", "noisbias": "noise bias of spectrum", } diff --git a/heracles/maps.py b/heracles/maps.py index 47439d2..cff8f16 100644 --- a/heracles/maps.py +++ b/heracles/maps.py @@ -331,9 +331,11 @@ def mapper(page: "CatalogPage") -> None: # set metadata of array update_metadata( pos, + catalog=catalog.label, spin=0, nbar=nbar, kernel="healpix", + nside=self.nside, power=power, bias=bias, ) @@ -434,9 +436,11 @@ def mapper(page: "CatalogPage") -> None: # set metadata of array update_metadata( val, + catalog=catalog.label, spin=0, wbar=wbar, kernel="healpix", + nside=self.nside, power=power, bias=bias, ) @@ -588,9 +592,11 @@ def mapper(page: "CatalogPage") -> None: # set metadata of array update_metadata( val, + catalog=catalog.label, spin=self.spin, wbar=wbar, kernel="healpix", + nside=self.nside, power=power, bias=bias, ) @@ -627,7 +633,14 @@ def __call__(self, catalog: "Catalog") -> MapData: # make a copy for updates to metadata vmap = np.copy(vmap) - update_metadata(vmap, spin=0, kernel="healpix", power=0) + update_metadata( + vmap, + catalog=catalog.label, + spin=0, + kernel="healpix", + nside=self.nside, + power=0, + ) return vmap @@ -691,7 +704,15 @@ def mapper(page: "CatalogPage") -> None: power = 1 # set metadata of arrays - update_metadata(wht, spin=0, wbar=wbar, kernel="healpix", power=power) + update_metadata( + wht, + catalog=catalog.label, + spin=0, + wbar=wbar, + kernel="healpix", + nside=self.nside, + power=power, + ) # return the weight map return wht @@ -849,7 +870,6 @@ def transform_maps( # convert maps to alms, taking care of complex and spin-weighted maps for (k, i), m in maps.items(): - nside = hp.get_nside(m) if isinstance(lmax, Mapping): _lmax = lmax.get((k, i)) or lmax.get(k) else: @@ -877,8 +897,7 @@ def transform_maps( alms = {(f"{k}_E", i): alms[1], (f"{k}_B", i): alms[2]} for ki, alm in alms.items(): - if md: - update_metadata(alm, nside=nside, **md) + update_metadata(alm, **md) out[ki] = alm del m, alms, alm diff --git a/tests/test_io.py b/tests/test_io.py index 46f4611..3e128cf 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -33,7 +33,15 @@ def mock_cls(rng): import numpy as np cl = rng.random(101) - cl.dtype = np.dtype(cl.dtype, metadata={"nside_1": 32, "nside_2": 64}) + cl.dtype = np.dtype( + cl.dtype, + metadata={ + "catalog_1": "cat-a.fits", + "nside_1": 32, + "catalog_2": "cat-b.fits", + "nside_2": 64, + }, + ) return { ("P", "P", 0, 0): cl, @@ -191,9 +199,9 @@ def test_write_read_maps(rng, tmp_path): v = rng.random(npix) g = rng.random((2, npix)) - p.dtype = np.dtype(p.dtype, metadata={"spin": 0}) - v.dtype = np.dtype(v.dtype, metadata={"spin": 0}) - g.dtype = np.dtype(g.dtype, metadata={"spin": 0}) + p.dtype = np.dtype(p.dtype, metadata={"catalog": "cat.fits", "spin": 0}) + v.dtype = np.dtype(v.dtype, metadata={"catalog": "cat.fits", "spin": 0}) + g.dtype = np.dtype(g.dtype, metadata={"catalog": "cat.fits", "spin": 2}) maps = { ("P", 1): p, diff --git a/tests/test_maps.py b/tests/test_maps.py index 60495d5..bfad8a5 100644 --- a/tests/test_maps.py +++ b/tests/test_maps.py @@ -100,7 +100,13 @@ def test_visibility_map(nside, vmap): assert result is not vmap assert result.shape == (12 * nside_out**2,) - assert result.dtype.metadata == {"spin": 0, "kernel": "healpix", "power": 0} + assert result.dtype.metadata == { + "catalog": catalog.label, + "spin": 0, + "kernel": "healpix", + "nside": nside_out, + "power": 0, + } assert np.isclose(result.mean(), fsky) # test missing visibility map @@ -124,15 +130,15 @@ def test_position_map(nside, catalog, vmap): nbar = 4.0 assert m.shape == (npix,) - assert m.dtype.metadata == pytest.approx( - { - "spin": 0, - "nbar": nbar, - "kernel": "healpix", - "power": 0, - "bias": bias / nbar**2, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 0, + "nbar": nbar, + "kernel": "healpix", + "nside": nside, + "power": 0, + "bias": pytest.approx(bias / nbar**2), + } np.testing.assert_array_equal(m, 0) # compute number count map @@ -140,15 +146,15 @@ def test_position_map(nside, catalog, vmap): m = map_catalog(PositionMap(nside, "ra", "dec", overdensity=False), catalog) assert m.shape == (npix,) - assert m.dtype.metadata == pytest.approx( - { - "spin": 0, - "nbar": 4.0, - "kernel": "healpix", - "power": 1, - "bias": bias, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 0, + "nbar": 4.0, + "kernel": "healpix", + "nside": nside, + "power": 1, + "bias": pytest.approx(bias), + } np.testing.assert_array_equal(m, 4) # compute overdensity maps with visibility map @@ -159,30 +165,30 @@ def test_position_map(nside, catalog, vmap): m = map_catalog(PositionMap(nside, "ra", "dec"), catalog) assert m.shape == (12 * nside**2,) - assert m.dtype.metadata == pytest.approx( - { - "spin": 0, - "nbar": nbar, - "kernel": "healpix", - "power": 0, - "bias": bias / nbar**2, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 0, + "nbar": pytest.approx(nbar), + "kernel": "healpix", + "nside": nside, + "power": 0, + "bias": pytest.approx(bias / nbar**2), + } # compute number count map with visibility map m = map_catalog(PositionMap(nside, "ra", "dec", overdensity=False), catalog) assert m.shape == (12 * nside**2,) - assert m.dtype.metadata == pytest.approx( - { - "spin": 0, - "nbar": nbar, - "kernel": "healpix", - "power": 1, - "bias": bias, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 0, + "nbar": pytest.approx(nbar), + "kernel": "healpix", + "nside": nside, + "power": 1, + "bias": pytest.approx(bias), + } def test_scalar_map(nside, catalog): @@ -200,29 +206,29 @@ def test_scalar_map(nside, catalog): bias = (4 * np.pi / npix / npix) * v2 assert m.shape == (npix,) - assert m.dtype.metadata == pytest.approx( - { - "spin": 0, - "wbar": wbar, - "kernel": "healpix", - "power": 0, - "bias": bias / wbar**2, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 0, + "wbar": pytest.approx(wbar), + "kernel": "healpix", + "nside": nside, + "power": 0, + "bias": pytest.approx(bias / wbar**2), + } np.testing.assert_array_almost_equal(m, 0) m = map_catalog(ScalarMap(nside, "ra", "dec", "g1", "w", normalize=False), catalog) assert m.shape == (npix,) - assert m.dtype.metadata == pytest.approx( - { - "spin": 0, - "wbar": wbar, - "kernel": "healpix", - "power": 1, - "bias": bias, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 0, + "wbar": pytest.approx(wbar), + "kernel": "healpix", + "nside": nside, + "power": 1, + "bias": pytest.approx(bias), + } np.testing.assert_array_almost_equal(m, 0) @@ -242,15 +248,15 @@ def test_complex_map(nside, catalog): bias = (4 * np.pi / npix / npix) * v2 / 2 assert m.shape == (2, npix) - assert m.dtype.metadata == pytest.approx( - { - "spin": 2, - "wbar": wbar, - "kernel": "healpix", - "power": 0, - "bias": bias / wbar**2, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 2, + "wbar": pytest.approx(wbar), + "kernel": "healpix", + "nside": nside, + "power": 0, + "bias": pytest.approx(bias / wbar**2), + } np.testing.assert_array_almost_equal(m, 0) m = map_catalog( @@ -259,15 +265,15 @@ def test_complex_map(nside, catalog): ) assert m.shape == (2, npix) - assert m.dtype.metadata == pytest.approx( - { - "spin": 1, - "wbar": wbar, - "kernel": "healpix", - "power": 1, - "bias": bias, - }, - ) + assert m.dtype.metadata == { + "catalog": catalog.label, + "spin": 1, + "wbar": pytest.approx(wbar), + "kernel": "healpix", + "nside": nside, + "power": 1, + "bias": pytest.approx(bias), + } np.testing.assert_array_almost_equal(m, 0) @@ -282,9 +288,11 @@ def test_weight_map(nside, catalog): assert m.shape == (12 * nside**2,) assert m.dtype.metadata == { + "catalog": catalog.label, "spin": 0, "wbar": wbar, "kernel": "healpix", + "nside": nside, "power": 0, } np.testing.assert_array_almost_equal(m, w / wbar) @@ -293,9 +301,11 @@ def test_weight_map(nside, catalog): assert m.shape == (12 * nside**2,) assert m.dtype.metadata == { + "catalog": catalog.label, "spin": 0, "wbar": wbar, "kernel": "healpix", + "nside": nside, "power": 1, } np.testing.assert_array_almost_equal(m, w) @@ -308,9 +318,9 @@ def test_transform_maps(rng): npix = 12 * nside**2 t = rng.standard_normal(npix) - update_metadata(t, spin=0, a=1) + update_metadata(t, spin=0, nside=nside, a=1) p = rng.standard_normal((2, npix)) - update_metadata(p, spin=2, b=2) + update_metadata(p, spin=2, nside=nside, b=2) # single scalar map maps = {("T", 0): t}