Skip to content

gh-216: move bias to spectra #417

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
JaimeRZP merged 3 commits into from
Jul 1, 2026
Merged

gh-216: move bias to spectra #417

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a26ab84
move bias to spectra
JaimeRZP Jul 1, 2026
2daffff
ruff
JaimeRZP Jul 1, 2026
3808daa
bug
JaimeRZP Jul 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 8 additions & 20 deletions heracles/fields.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,15 +304,12 @@ 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
)
update_metadata(pos, catalog, nbar=nbar, musq=musq, dens=dens, fsky=fsky)

# return the position map
return pos
Expand Down Expand Up @@ -375,16 +372,13 @@ 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
)
update_metadata(val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky)

# return the value map
return val
Expand Down Expand Up @@ -451,16 +445,13 @@ 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
)
update_metadata(val, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky)

# return the shear map
return val
Expand Down Expand Up @@ -556,16 +547,13 @@ 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
)
update_metadata(wht, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky)

# return the weight map
return wht
Expand Down
21 changes: 11 additions & 10 deletions heracles/twopoint.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Comment on lines +259 to +267
if bias is not None:
md["bias"] = bias

Expand Down
15 changes: 0 additions & 15 deletions tests/test_fields.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand All @@ -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)

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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):
Expand All @@ -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,
Expand All @@ -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)

Expand All @@ -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,
Expand All @@ -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)

Expand All @@ -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)
Expand All @@ -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)

Expand Down
42 changes: 42 additions & 0 deletions tests/test_twopoint.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,48 @@ 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
)
Comment on lines +166 to +181

# cross-spectrum: no bias key
alms_cross = {("F", 0): a, ("F", 1): a.copy()}
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
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

Expand Down
Loading