Skip to content

Commit

Permalink
ENH: implement theoretical bias values (#22)
Browse files Browse the repository at this point in the history 
This change implements theoretical values for the additive biases, which
are reported in the `"bias"` metadata keyword. Also removes usage of
"noise bias" as the generic term, since only some of the additive biases
are of a stochastic nature.
  • Loading branch information
@ntessore
ntessore authored Sep 6, 2023
1 parent 4f41666 commit 624dcbe
Show file tree
Hide file tree
Showing 5 changed files with 255 additions and 216 deletions.
181 changes: 51 additions & 130 deletions examples/example.ipynb
View file

Large diffs are not rendered by default.

98 changes: 85 additions & 13 deletions heracles/maps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,9 +266,11 @@ def mapper(page: "CatalogPage") -> None:

ngal += page.size

# call yields here to apply mapper over entire catalogue
# the mapper function is yield-ed to be applied over the catalogue
yield mapper

# when function resumes, mapping has finished

# get visibility map if present in catalogue
vmap = catalog.visibility

Expand All @@ -287,8 +289,11 @@ def mapper(page: "CatalogPage") -> None:

# compute average number density
nbar = ngal / npix
if vmap is not None:
nbar /= np.mean(vmap)
if vmap is None:
vbar = 1
else:
vbar = np.mean(vmap)
nbar /= vbar

# compute overdensity if asked to
if self._overdensity:
Expand All @@ -298,11 +303,20 @@ def mapper(page: "CatalogPage") -> None:
else:
pos -= vmap
power = 0
bias = 4 * np.pi * vbar**2 / ngal
else:
power = 1
bias = (4 * np.pi / npix) * (ngal / npix)

# set metadata of array
update_metadata(pos, spin=0, nbar=nbar, kernel="healpix", power=power)
update_metadata(
pos,
spin=0,
nbar=nbar,
kernel="healpix",
power=power,
bias=bias,
)

# return the position map
return pos
Expand Down Expand Up @@ -343,8 +357,14 @@ def __call__(self, catalog: "Catalog") -> MapGenerator:
wht = np.zeros(npix)
val = np.zeros(npix)

# total weighted variance from online algorithm
ngal = 0
wmean, var = 0.0, 0.0

# go through pages in catalogue and map values
def mapper(page: "CatalogPage") -> None:
nonlocal ngal, wmean, var

if wcol is not None:
page.delete(page[wcol] == 0)

Expand All @@ -359,25 +379,46 @@ def mapper(page: "CatalogPage") -> None:

_map_real(wht, val, ipix, w, v)

# call yields here to apply mapper over entire catalogue
ngal += page.size
wmean += (w - wmean).sum() / ngal
var += ((w * v) ** 2 - var).sum() / ngal

# the mapper function is yield-ed to be applied over the catalogue
yield mapper

# compute average weight in nonzero pixels
wbar = wht.mean()
# when function resumes, mapping has finished

# compute mean visibility
if catalog.visibility is None:
vbar = 1
else:
vbar = np.mean(catalog.visibility)

# compute mean weight per visible pixel
wbar = ngal / npix / vbar * wmean

# normalise the weight in each pixel if asked to
if self.normalize:
wht /= wbar
power = 0
bias = 4 * np.pi * vbar**2 / ngal * (var / wmean**2)
else:
power = 1
bias = (4 * np.pi / npix) * (ngal / npix) * var

# value was averaged in each pixel for numerical stability
# now compute the sum
val *= wht

# set metadata of array
update_metadata(val, spin=0, wbar=wbar, kernel="healpix", power=power)
update_metadata(
val,
spin=0,
wbar=wbar,
kernel="healpix",
power=power,
bias=bias,
)

# return the value map
return val
Expand Down Expand Up @@ -457,9 +498,15 @@ def __call__(self, catalog: "Catalog") -> MapGenerator:
wht = np.zeros(npix)
val = np.zeros((2, npix))

# total weighted variance from online algorithm
ngal = 0
wmean, var = 0.0, 0.0

# go through pages in catalogue and get the shear values,
# randomise if asked to, and do the mapping
def mapper(page: "CatalogPage") -> None:
nonlocal ngal, wmean, var

if wcol is not None:
page.delete(page[wcol] == 0)

Expand All @@ -483,25 +530,46 @@ def mapper(page: "CatalogPage") -> None:

_map_complex(wht, val, ipix, w, re, im)

# call yields here to apply mapper over entire catalogue
ngal += page.size
wmean += (w - wmean).sum() / ngal
var += ((w * re) ** 2 + (w * im) ** 2 - var).sum() / ngal

# the mapper function is yield-ed to be applied over the catalogue
yield mapper

# compute average weight in nonzero pixels
wbar = wht.mean()
# when function resumes, mapping has finished

# compute mean visibility
if catalog.visibility is None:
vbar = 1
else:
vbar = np.mean(catalog.visibility)

# mean weight per visible pixel
wbar = ngal / npix / vbar * wmean

# normalise the weight in each pixel if asked to
if self.normalize:
wht /= wbar
power = 0
bias = 2 * np.pi * vbar**2 / ngal * (var / wmean**2)
else:
power = 1
bias = (2 * np.pi / npix) * (ngal / npix) * var

# value was averaged in each pixel for numerical stability
# now compute the sum
val *= wht

# set metadata of array
update_metadata(val, spin=self.spin, wbar=wbar, kernel="healpix", power=power)
update_metadata(
val,
spin=self.spin,
wbar=wbar,
kernel="healpix",
power=power,
bias=bias,
)

# return the shear map
return val
Expand Down Expand Up @@ -581,11 +649,15 @@ def mapper(page: "CatalogPage") -> None:

_map_weight(wht, ipix, w)

# call yields here to apply mapper over entire catalogue
# the mapper function is yield-ed to be applied over the catalogue
yield mapper

# when function resumes, mapping has finished

# compute average weight in nonzero pixels
wbar = wht.mean()
if catalog.visibility is not None:
wbar /= np.mean(catalog.visibility)

# normalise the weight in each pixel if asked to
if self.normalize:
Expand Down
38 changes: 19 additions & 19 deletions heracles/twopoint.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,13 +88,13 @@ def angular_power_spectra(alms, alms2=None, *, lmax=None, include=None, exclude=
md = {}
if alm1.dtype.metadata:
for key, value in alm1.dtype.metadata.items():
if key == "noisbias":
if key == "bias":
md[key] = value if k1 == k2 and i1 == i2 else 0.0
else:
md[f"{key}_1"] = value
if alm2.dtype.metadata:
for key, value in alm2.dtype.metadata.items():
if key == "noisbias":
if key == "bias":
pass
else:
md[f"{key}_2"] = value
Expand All @@ -116,19 +116,16 @@ def angular_power_spectra(alms, alms2=None, *, lmax=None, include=None, exclude=
return cls


def debias_cls(cls, noisebias=None, *, inplace=False):
"""remove noise bias from cls"""
def debias_cls(cls, bias=None, *, inplace=False):
"""remove bias from cls"""

logger.info("debiasing %d cl(s)%s", len(cls), " in place" if inplace else "")
t = time.monotonic()

# toc dict for noise biases
nbs = noisebias or {}

# the output toc dict
out = cls if inplace else {}

# subtract noise bias of each cl in turn
# subtract bias of each cl in turn
for key in cls:
logger.info("debiasing %s cl for bins %s, %s", *key)

Expand All @@ -142,14 +139,17 @@ def debias_cls(cls, noisebias=None, *, inplace=False):
# minimum l for correction
lmin = max(abs(md.get("spin_1", 0)), abs(md.get("spin_2", 0)))

# get noise bias from explicit dict, if given, or metadata
nb = nbs.get(key, md.get("noisbias", 0.0))
# get bias from explicit dict, if given, or metadata
if bias is None:
b = md.get("bias", 0.0)
else:
b = bias.get(key, 0.0)

# remove noise bias
cl[lmin:] -= nb
# remove bias
cl[lmin:] -= b

# write noise bias to corrected cl
update_metadata(cl, noisbias=nb)
update_metadata(cl, bias=b)

# store debiased cl in output set
out[key] = cl
Expand Down Expand Up @@ -425,7 +425,7 @@ def norm(a, b):
return out


def random_noisebias(
def random_bias(
maps,
catalogs,
*,
Expand All @@ -437,13 +437,13 @@ def random_noisebias(
progress=False,
**kwargs,
):
"""noise bias estimate from randomised position and shear maps
"""bias estimate from randomised maps
The ``include`` and ``exclude`` selection is applied to the maps.
"""

logger.info("estimating two-point noise bias for %d catalog(s)", len(catalogs))
logger.info("estimating two-point bias for %d catalog(s)", len(catalogs))
logger.info("randomising %s maps", ", ".join(map(str, maps)))
t = time.monotonic()

Expand All @@ -453,7 +453,7 @@ def random_noisebias(
# include will be set below after we have the first set of alms
include_cls = None
if full:
logger.info("estimating cross-noise biases")
logger.info("estimating cross-biases")

# set all input maps to randomize
# store and later reset their initial state
Expand All @@ -466,7 +466,7 @@ def random_noisebias(

for n in range(repeat):
logger.info(
"estimating noise bias from randomised maps%s",
"estimating bias from randomised maps%s",
"" if n == 0 else f" (repeat {n+1})",
)

Expand Down Expand Up @@ -495,7 +495,7 @@ def random_noisebias(
m.randomize = randomize[k]

logger.info(
"estimated %d two-point noise biases in %s",
"estimated %d two-point biases in %s",
len(nbs),
timedelta(seconds=(time.monotonic() - t)),
)
Expand Down
Loading

0 comments on commit 624dcbe

Please sign in to comment.