From a26ab84b87743e63a684d3a5766b151ba2b3d75c Mon Sep 17 00:00:00 2001 From: jaimerzp Date: Wed, 1 Jul 2026 11:50:01 +0100 Subject: [PATCH 1/3] move bias to spectra --- heracles/fields.py | 20 ++++++++------------ heracles/twopoint.py | 21 +++++++++++---------- tests/test_fields.py | 15 --------------- tests/test_twopoint.py | 34 ++++++++++++++++++++++++++++++++++ 4 files changed, 53 insertions(+), 37 deletions(-) diff --git a/heracles/fields.py b/heracles/fields.py index 31d5ea2a..d4fea4c7 100644 --- a/heracles/fields.py +++ b/heracles/fields.py @@ -304,14 +304,13 @@ async def __call__( pos -= vis del vis - # compute bias of positions, including weight variance + # compute bias ingredients for positions, including weight variance musq = 1.0 dens = (nbar / mapper.area) ** 2 / (ngal / (4 * np.pi * fsky)) / w2mean - bias = fsky * musq / dens # set metadata of array update_metadata( - pos, catalog, nbar=nbar, musq=musq, dens=dens, fsky=fsky, bias=bias + pos, catalog, nbar=nbar, musq=musq, dens=dens, fsky=fsky ) # return the position map @@ -375,15 +374,14 @@ async def __call__( # normalise the map val /= wbar - # compute bias from variance (per object) + # compute bias ingredients from variance (per object) musq = var / w2mean deff = w2mean / wmean**2 dens = ngal / (4 * np.pi * fsky) / deff - bias = fsky * musq / dens # set metadata of array update_metadata( - val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky, bias=bias + val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky ) # return the value map @@ -451,15 +449,14 @@ async def __call__( # normalise the map val /= wbar - # bias from measured variance, for E/B decomposition + # bias ingredients from measured variance, for E/B decomposition musq = var / w2mean deff = w2mean / wmean**2 dens = ngal / (4 * np.pi * fsky) / deff - bias = (1 / 2) * fsky * musq / dens # set metadata of array update_metadata( - val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky, bias=bias + val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky ) # return the shear map @@ -556,15 +553,14 @@ async def __call__( # normalise the map wht /= wbar - # bias from weights + # bias ingredients from weights musq = 1.0 deff = w2mean / wmean**2 dens = ngal / (4 * np.pi * fsky) / deff - bias = fsky * musq / dens # set metadata of array update_metadata( - wht, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky, bias=bias + wht, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky ) # return the weight map diff --git a/heracles/twopoint.py b/heracles/twopoint.py index 1472095e..0d1a0b1c 100644 --- a/heracles/twopoint.py +++ b/heracles/twopoint.py @@ -252,18 +252,19 @@ def angular_power_spectra( raise ValueError(f"missing spin metadata for {k1} or {k2}") # collect metadata md = {} - bias = None for key, value in md1.items(): - if key == "bias": - if k1 == k2 and i1 == i2: - bias = value - else: - md[f"{key}_1"] = value + md[f"{key}_1"] = value for key, value in md2.items(): - if key == "bias": - pass - else: - md[f"{key}_2"] = value + md[f"{key}_2"] = value + # compute bias for auto-spectra from ingredients stored during field mapping + bias = None + if k1 == k2 and i1 == i2: + _fsky = md1.get("fsky") + _musq = md1.get("musq") + _dens = md1.get("dens") + if _fsky is not None and _musq is not None and _dens is not None: + factor = 0.5 if s1 == s2 == 2 else 1.0 + bias = factor * _fsky * _musq / _dens if bias is not None: md["bias"] = bias diff --git a/tests/test_fields.py b/tests/test_fields.py index e4909b8a..e9de34e1 100644 --- a/tests/test_fields.py +++ b/tests/test_fields.py @@ -172,9 +172,7 @@ def test_visibility(nside, vmap): def test_positions(mapper, catalog, vmap): from heracles.fields import Positions - # bias npix = 12 * mapper.nside**2 - bias = (4 * np.pi / npix) * (catalog.size / npix) # normal mode: compute overdensity maps with metadata @@ -207,7 +205,6 @@ def test_positions(mapper, catalog, vmap): "musq": 1.0, "dens": pytest.approx(npix / np.pi), "fsky": 1.0, - "bias": pytest.approx(bias / nbar**2), } np.testing.assert_array_equal(m, 0) @@ -229,7 +226,6 @@ def test_positions(mapper, catalog, vmap): "musq": 1.0, "dens": pytest.approx(npix / np.pi), "fsky": 1.0, - "bias": pytest.approx(bias / nbar**2), } np.testing.assert_array_equal(m, 1.0) @@ -255,7 +251,6 @@ def test_positions(mapper, catalog, vmap): "musq": 1.0, "dens": pytest.approx(npix / (np.pi * catalog.fsky)), "fsky": catalog.fsky, - "bias": pytest.approx(bias / nbar**2), } # compute number count map with visibility map @@ -276,7 +271,6 @@ def test_positions(mapper, catalog, vmap): "musq": 1.0, "dens": pytest.approx(npix / (np.pi * catalog.fsky)), "fsky": catalog.fsky, - "bias": pytest.approx(bias / nbar**2), } # compute overdensity maps with given (incorrect) nbar @@ -286,7 +280,6 @@ def test_positions(mapper, catalog, vmap): m = coroutines.run(f(catalog)) assert m.dtype.metadata["nbar"] == 2 * nbar - assert m.dtype.metadata["bias"] == pytest.approx(bias / (2 * nbar) ** 2) def test_scalar_field(mapper, catalog): @@ -310,8 +303,6 @@ def test_scalar_field(mapper, catalog): musq = var / w2mean deff = w2mean / wmean**2 dens = npix / np.pi / deff - bias = (4 * np.pi / npix / npix) * v2 - assert m.shape == (npix,) assert m.dtype.metadata == { "catalog": catalog.label, @@ -325,7 +316,6 @@ def test_scalar_field(mapper, catalog): "musq": pytest.approx(musq), "dens": pytest.approx(dens), "fsky": 1.0, - "bias": pytest.approx(bias / wbar**2), } np.testing.assert_array_almost_equal(m, 0) @@ -352,8 +342,6 @@ def test_complex_field(mapper, catalog): musq = var / w2mean deff = w2mean / wmean**2 dens = npix / np.pi / deff - bias = (4 * np.pi / npix / npix) * v2 / 2 - assert m.shape == (2, npix) assert m.dtype.metadata == { "catalog": catalog.label, @@ -367,7 +355,6 @@ def test_complex_field(mapper, catalog): "fsky": 1.0, "dens": pytest.approx(dens), "musq": pytest.approx(musq), - "bias": pytest.approx(bias / wbar**2), } np.testing.assert_array_almost_equal(m, 0) @@ -384,7 +371,6 @@ def test_weights(mapper, catalog): v2 = (w**2).sum() w = w.reshape(w.size // 4, 4).sum(axis=-1) wbar = w.mean() - bias = (4 * np.pi / npix / npix) * v2 v1 = w.sum() wmean = v1 / (4.0 * npix) w2mean = v2 / (4.0 * npix) @@ -404,7 +390,6 @@ def test_weights(mapper, catalog): "musq": 1.0, "dens": pytest.approx(dens), "fsky": 1.0, - "bias": pytest.approx(bias / wbar**2), } np.testing.assert_array_almost_equal(m, w / wbar) diff --git a/tests/test_twopoint.py b/tests/test_twopoint.py index 50d59bce..740e42be 100644 --- a/tests/test_twopoint.py +++ b/tests/test_twopoint.py @@ -151,6 +151,40 @@ def test_angular_power_spectra(mock_alms, lmax): assert next(call_iter) == call((a, b, i, j), inc, exc) +def test_angular_power_spectra_bias(lmax): + from heracles.twopoint import angular_power_spectra + + size = (lmax + 1) * (lmax + 2) // 2 + + fsky, musq, dens = 0.5, 1.2, 3.4 + + # spin-0 auto: bias = fsky * musq / dens + a = np.zeros(size, dtype=complex) + a.dtype = np.dtype(a.dtype, metadata={"spin": 0, "fsky": fsky, "musq": musq, "dens": dens}) + alms = {("F", 0): a} + cls = angular_power_spectra(alms, debias=False) + assert cls["F", "F", 0, 0].dtype.metadata["bias"] == pytest.approx(fsky * musq / dens) + + # spin-2 auto: bias = 0.5 * fsky * musq / dens + b = np.zeros((2, size), dtype=complex) + b.dtype = np.dtype(b.dtype, metadata={"spin": 2, "fsky": fsky, "musq": musq, "dens": dens}) + alms2 = {("G", 0): b} + cls2 = angular_power_spectra(alms2, debias=False) + assert cls2["G", "G", 0, 0].dtype.metadata["bias"] == pytest.approx(0.5 * fsky * musq / dens) + + # cross-spectrum: no bias key + alms_cross = {("F", 0): a, ("F", 1): a.copy()} + cls_cross = angular_power_spectra({("F", 0): a}, {("F", 1): a.copy()}, debias=False) + assert "bias" not in (cls_cross["F", "F", 0, 1].dtype.metadata or {}) + + # external map (no ingredients): no bias key + c = np.zeros(size, dtype=complex) + c.dtype = np.dtype(c.dtype, metadata={"spin": 0}) + alms_ext = {("F", 0): c} + cls_ext = angular_power_spectra(alms_ext, debias=False) + assert "bias" not in (cls_ext["F", "F", 0, 0].dtype.metadata or {}) + + def test_debias_cls(): from heracles.twopoint import debias_cls From 2daffff299e6778595c0c22fdb19f216f110d998 Mon Sep 17 00:00:00 2001 From: jaimerzp Date: Wed, 1 Jul 2026 11:52:19 +0100 Subject: [PATCH 2/3] ruff --- heracles/fields.py | 16 ++++------------ tests/test_twopoint.py | 16 ++++++++++++---- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/heracles/fields.py b/heracles/fields.py index d4fea4c7..aa68ad12 100644 --- a/heracles/fields.py +++ b/heracles/fields.py @@ -309,9 +309,7 @@ async def __call__( dens = (nbar / mapper.area) ** 2 / (ngal / (4 * np.pi * fsky)) / w2mean # set metadata of array - update_metadata( - pos, catalog, nbar=nbar, musq=musq, dens=dens, fsky=fsky - ) + update_metadata(pos, catalog, nbar=nbar, musq=musq, dens=dens, fsky=fsky) # return the position map return pos @@ -380,9 +378,7 @@ async def __call__( dens = ngal / (4 * np.pi * fsky) / deff # set metadata of array - update_metadata( - val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky - ) + update_metadata(val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky) # return the value map return val @@ -455,9 +451,7 @@ async def __call__( dens = ngal / (4 * np.pi * fsky) / deff # set metadata of array - update_metadata( - val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky - ) + update_metadata(val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky) # return the shear map return val @@ -559,9 +553,7 @@ async def __call__( dens = ngal / (4 * np.pi * fsky) / deff # set metadata of array - update_metadata( - wht, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky - ) + update_metadata(wht, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky) # return the weight map return wht diff --git a/tests/test_twopoint.py b/tests/test_twopoint.py index 740e42be..1fd10df5 100644 --- a/tests/test_twopoint.py +++ b/tests/test_twopoint.py @@ -160,17 +160,25 @@ def test_angular_power_spectra_bias(lmax): # spin-0 auto: bias = fsky * musq / dens a = np.zeros(size, dtype=complex) - a.dtype = np.dtype(a.dtype, metadata={"spin": 0, "fsky": fsky, "musq": musq, "dens": dens}) + a.dtype = np.dtype( + a.dtype, metadata={"spin": 0, "fsky": fsky, "musq": musq, "dens": dens} + ) alms = {("F", 0): a} cls = angular_power_spectra(alms, debias=False) - assert cls["F", "F", 0, 0].dtype.metadata["bias"] == pytest.approx(fsky * musq / dens) + assert cls["F", "F", 0, 0].dtype.metadata["bias"] == pytest.approx( + fsky * musq / dens + ) # spin-2 auto: bias = 0.5 * fsky * musq / dens b = np.zeros((2, size), dtype=complex) - b.dtype = np.dtype(b.dtype, metadata={"spin": 2, "fsky": fsky, "musq": musq, "dens": dens}) + b.dtype = np.dtype( + b.dtype, metadata={"spin": 2, "fsky": fsky, "musq": musq, "dens": dens} + ) alms2 = {("G", 0): b} cls2 = angular_power_spectra(alms2, debias=False) - assert cls2["G", "G", 0, 0].dtype.metadata["bias"] == pytest.approx(0.5 * fsky * musq / dens) + assert cls2["G", "G", 0, 0].dtype.metadata["bias"] == pytest.approx( + 0.5 * fsky * musq / dens + ) # cross-spectrum: no bias key alms_cross = {("F", 0): a, ("F", 1): a.copy()} From 3808daa610d083c26dc09ade5591b5faafeaf6df Mon Sep 17 00:00:00 2001 From: jaimerzp Date: Wed, 1 Jul 2026 11:54:02 +0100 Subject: [PATCH 3/3] bug --- tests/test_twopoint.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_twopoint.py b/tests/test_twopoint.py index 1fd10df5..3a050be5 100644 --- a/tests/test_twopoint.py +++ b/tests/test_twopoint.py @@ -182,7 +182,7 @@ def test_angular_power_spectra_bias(lmax): # cross-spectrum: no bias key alms_cross = {("F", 0): a, ("F", 1): a.copy()} - cls_cross = angular_power_spectra({("F", 0): a}, {("F", 1): a.copy()}, debias=False) + cls_cross = angular_power_spectra(alms_cross, debias=False) assert "bias" not in (cls_cross["F", "F", 0, 1].dtype.metadata or {}) # external map (no ingredients): no bias key