From 06ffc27ad61478c68fe869868b2d5b4b399b81f1 Mon Sep 17 00:00:00 2001 From: Nicolas Tessore Date: Sat, 6 Jan 2024 10:07:59 +0000 Subject: [PATCH] ENH(twopoint): debias spectra at time of computation (#104) 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 --- examples/example.ipynb | 108 +++++++++++++++-------------------------- heracles/twopoint.py | 7 +++ 2 files changed, 45 insertions(+), 70 deletions(-) diff --git a/examples/example.ipynb b/examples/example.ipynb index a70da86..3cbc5b0 100644 --- a/examples/example.ipynb +++ b/examples/example.ipynb @@ -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." ] }, { @@ -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", "}" ] }, @@ -522,7 +522,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "69f552bc18b543a5a39235bda3703ff2", + "model_id": "8dd05200ab024231be455bd2d581fa95", "version_major": 2, "version_minor": 0 }, @@ -536,13 +536,13 @@ { "data": { "text/html": [ - "
/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",
        "
\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"
       ]
@@ -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."
    ]
   },
   {
@@ -692,7 +692,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "a78c0a381db24cfda4bf1521bd78fc4f",
+       "model_id": "bd982490753c4c03bcacaf0551fba9d7",
        "version_major": 2,
        "version_minor": 0
       },
@@ -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)."
    ]
   },
   {
@@ -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": {},
@@ -893,7 +861,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 27,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -910,7 +878,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -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
       },
@@ -982,7 +950,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 30,
    "metadata": {},
    "outputs": [
     {
@@ -1005,7 +973,7 @@
        " ('V', 7)]"
       ]
      },
-     "execution_count": 32,
+     "execution_count": 30,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -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
       },
@@ -1077,7 +1045,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 32,
    "metadata": {
     "scrolled": true
    },
@@ -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": [
@@ -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
       },
@@ -1173,7 +1141,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 35,
    "metadata": {},
    "outputs": [
     {
@@ -1194,7 +1162,7 @@
        " ('G_E', 'G_B', 9, 9)]"
       ]
      },
-     "execution_count": 37,
+     "execution_count": 35,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1219,7 +1187,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 36,
    "metadata": {},
    "outputs": [
     {
@@ -1256,7 +1224,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 37,
    "metadata": {},
    "outputs": [
     {
@@ -1293,7 +1261,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 40,
+   "execution_count": 38,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1319,7 +1287,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 41,
+   "execution_count": 39,
    "metadata": {},
    "outputs": [
     {
@@ -1478,7 +1446,7 @@
        " "
       ]
      },
-     "execution_count": 41,
+     "execution_count": 39,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1501,7 +1469,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 42,
+   "execution_count": 40,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1517,7 +1485,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
+   "execution_count": 41,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1534,7 +1502,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 44,
+   "execution_count": 42,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1557,7 +1525,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 45,
+   "execution_count": 43,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1574,7 +1542,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 46,
+   "execution_count": 44,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1592,7 +1560,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 47,
+   "execution_count": 45,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1634,7 +1602,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 48,
+   "execution_count": 46,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1643,7 +1611,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 49,
+   "execution_count": 47,
    "metadata": {},
    "outputs": [
     {
@@ -1668,7 +1636,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 50,
+   "execution_count": 48,
    "metadata": {},
    "outputs": [
     {
diff --git a/heracles/twopoint.py b/heracles/twopoint.py
index e630026..2857801 100644
--- a/heracles/twopoint.py
+++ b/heracles/twopoint.py
@@ -51,6 +51,7 @@ def angular_power_spectra(
     alms2=None,
     *,
     lmax=None,
+    debias=True,
     bins=None,
     weights=None,
     include=None,
@@ -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)