Skip to content

Commit

Permalink
ENH(twopoint): debias spectra at time of computation (#104)
Browse files Browse the repository at this point in the history 
Add a `debias=True` parameter to `angular_power_spectra()` which causes
the noise bias to be removed from spectra as soon as they are computed.

Closes: #24
  • Loading branch information
@ntessore
ntessore authored Jan 6, 2024
1 parent 115e1eb commit 06ffc27
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 70 deletions.
108 changes: 38 additions & 70 deletions examples/example.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -454,7 +454,7 @@
"\n",
"We also specify that the shear field should flip the sign of the \"G2\" column by adding a minus sign.\n",
"\n",
"Finally, we define the optional names of the weight functions (`\"V\"`, `\"W\"`) of the fields. These will be used further down below to compute mixing matrices for the fields."
"Finally, we define the optional names of the masks (`\"V\"`, `\"W\"`) of the fields. These will be used further down below to compute mixing matrices."
]
},
{
Expand All @@ -467,8 +467,8 @@
"lonlat = ('RIGHT_ASCENSION', 'DECLINATION')\n",
"\n",
"fields = {\n",
" 'P': Positions(*lonlat, weight=\"V\"),\n",
" 'G': Shears(*lonlat, 'G1', '-G2', 'WEIGHT', weight=\"W\"),\n",
" 'P': Positions(*lonlat, mask=\"V\"),\n",
" 'G': Shears(*lonlat, 'G1', '-G2', 'WEIGHT', mask=\"W\"),\n",
"}"
]
},
Expand Down Expand Up @@ -522,7 +522,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "69f552bc18b543a5a39235bda3703ff2",
"model_id": "8dd05200ab024231be455bd2d581fa95",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -536,13 +536,13 @@
{
"data": {
"text/html": [
"<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">/Users/ntessore/code/heracles-ec/heracles/heracles/fields.py:254: UserWarning: position and visibility maps have \n",
"<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">/Users/ntessore/code/heracles-ec/heracles/heracles/fields.py:256: UserWarning: position and visibility maps have \n",
"different NSIDE\n",
" warnings.warn(\"position and visibility maps have different NSIDE\")\n",
"</pre>\n"
],
"text/plain": [
"/Users/ntessore/code/heracles-ec/heracles/heracles/fields.py:254: UserWarning: position and visibility maps have \n",
"/Users/ntessore/code/heracles-ec/heracles/heracles/fields.py:256: UserWarning: position and visibility maps have \n",
"different NSIDE\n",
" warnings.warn(\"position and visibility maps have different NSIDE\")\n"
]
Expand Down Expand Up @@ -672,7 +672,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we are working in harmonic space, the real-space maps we created so far are not what we are actually after. We therefore transform them into harmonic-space maps (i.e. $a_{lm}$) using the `transform_maps()` function."
"Since we are working in harmonic space, the real-space maps we created so far are not what we are actually after. We therefore transform them into harmonic-space maps (i.e. $a_{lm}$) using the `transform_maps()` function. The function will also remove the HEALPix pixel window function, unless `deconvolve=False` is passed."
]
},
{
Expand All @@ -692,7 +692,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "a78c0a381db24cfda4bf1521bd78fc4f",
"model_id": "bd982490753c4c03bcacaf0551fba9d7",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -778,7 +778,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"We are now able to compute the two-point statistics in the form of angular power spectra. We simply call the `angular_power_spectrum()` function on the `alms`."
"We are now able to compute the two-point statistics in the form of angular power spectra. We simply call the `angular_power_spectrum()` function on the `alms`. The function will also remove the noise bias from the spectra (unless `debias=False` is passed), and can optionally compute binned spectra (using the `bins=` and `weights=` parameters)."
]
},
{
Expand Down Expand Up @@ -838,38 +838,6 @@
"list(cls.keys())[:6] + ['...'] + list(cls.keys())[-6:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Debiasing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can debias our angular power spectra by subtracting the bias terms from the auto-correlations. Since we do not need the biased `cls`, we compute the debiased ones in place."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"from heracles.twopoint import debias_cls"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"debias_cls(cls, inplace=True);"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand All @@ -893,7 +861,7 @@
},
{
"cell_type": "code",
"execution_count": 29,
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -910,7 +878,7 @@
},
{
"cell_type": "code",
"execution_count": 30,
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -919,13 +887,13 @@
},
{
"cell_type": "code",
"execution_count": 31,
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6e4d3dc39afc46a9b31d03ccd3fb7b14",
"model_id": "0345e186017a46b8916f43e440b12ee7",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -982,7 +950,7 @@
},
{
"cell_type": "code",
"execution_count": 32,
"execution_count": 30,
"metadata": {},
"outputs": [
{
Expand All @@ -1005,7 +973,7 @@
" ('V', 7)]"
]
},
"execution_count": 32,
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -1023,13 +991,13 @@
},
{
"cell_type": "code",
"execution_count": 33,
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "34f59966f570417880c30c9fd306d745",
"model_id": "f0a98718f7e24414bb9076c6805eed92",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -1077,7 +1045,7 @@
},
{
"cell_type": "code",
"execution_count": 34,
"execution_count": 32,
"metadata": {
"scrolled": true
},
Expand All @@ -1090,12 +1058,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"With the angular power spectra of visibility and weight maps available, we can compute all mixing matrices with the `mixing_matrices()` function."
"With the angular power spectra of visibility and weight maps available, we can compute all mixing matrices with the `mixing_matrices()` function. The function can optionally compute mixing matrices for binned spectra using the `bins=` and `weights=` parameters."
]
},
{
"cell_type": "code",
"execution_count": 35,
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1104,13 +1072,13 @@
},
{
"cell_type": "code",
"execution_count": 36,
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "230f11f3d8e741858a241369d682cc43",
"model_id": "6e7bc7d249504a888c6249089ff46b7f",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -1173,7 +1141,7 @@
},
{
"cell_type": "code",
"execution_count": 37,
"execution_count": 35,
"metadata": {},
"outputs": [
{
Expand All @@ -1194,7 +1162,7 @@
" ('G_E', 'G_B', 9, 9)]"
]
},
"execution_count": 37,
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -1219,7 +1187,7 @@
},
{
"cell_type": "code",
"execution_count": 38,
"execution_count": 36,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -1256,7 +1224,7 @@
},
{
"cell_type": "code",
"execution_count": 39,
"execution_count": 37,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -1293,7 +1261,7 @@
},
{
"cell_type": "code",
"execution_count": 40,
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1319,7 +1287,7 @@
},
{
"cell_type": "code",
"execution_count": 41,
"execution_count": 39,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -1478,7 +1446,7 @@
" "
]
},
"execution_count": 41,
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -1501,7 +1469,7 @@
},
{
"cell_type": "code",
"execution_count": 42,
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1517,7 +1485,7 @@
},
{
"cell_type": "code",
"execution_count": 43,
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1534,7 +1502,7 @@
},
{
"cell_type": "code",
"execution_count": 44,
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1557,7 +1525,7 @@
},
{
"cell_type": "code",
"execution_count": 45,
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1574,7 +1542,7 @@
},
{
"cell_type": "code",
"execution_count": 46,
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1592,7 +1560,7 @@
},
{
"cell_type": "code",
"execution_count": 47,
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -1634,7 +1602,7 @@
},
{
"cell_type": "code",
"execution_count": 48,
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -1643,7 +1611,7 @@
},
{
"cell_type": "code",
"execution_count": 49,
"execution_count": 47,
"metadata": {},
"outputs": [
{
Expand All @@ -1668,7 +1636,7 @@
},
{
"cell_type": "code",
"execution_count": 50,
"execution_count": 48,
"metadata": {},
"outputs": [
{
Expand Down
7 changes: 7 additions & 0 deletions heracles/twopoint.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ def angular_power_spectra(
alms2=None,
*,
lmax=None,
debias=True,
bins=None,
weights=None,
include=None,
Expand Down Expand Up @@ -131,6 +132,12 @@ def angular_power_spectra(
md[f"{key}_2"] = value
update_metadata(cl, **md)

# debias cl if asked to
if debias:
# minimum l for correction
_lmin = max(abs(md.get("spin_1", 0)), abs(md.get("spin_2", 0)))
cl[_lmin:] -= md.get("bias", 0.0)

# if bins are given, apply the binning
if bins is not None:
cl = bin2pt(cl, bins, "CL", weights=weights)
Expand Down

0 comments on commit 06ffc27

Please sign in to comment.