working weghts

heracles-ec · Nov 27, 2024 · 573d557 · 573d557
1 parent 2288f0a
commit 573d557
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 15 deletions.
diff --git a/heracles/fields.py b/heracles/fields.py
@@ -306,7 +306,7 @@ async def __call__(
 
         # compute bias of positions, including weight variance
         musq = 1.0
-        dens = (nbar / mapper.area)**2 / (ngal / (4 * np.pi * fsky)) / w2mean
+        dens = (nbar / mapper.area) ** 2 / (ngal / (4 * np.pi * fsky)) / w2mean
         bias = fsky * musq / dens
 
         # set metadata of array
@@ -377,9 +377,9 @@ async def __call__(
         # compute bias from variance (per object)
         musq = var / wmean**2
         dens = ngal / (4 * np.pi * fsky)
-        #variance = var / w2mean
-        #deff = w2mean / wmean**2
-        #neff = ngal / (4 * np.pi * fsky) / deff
+        # variance = var / w2mean
+        # deff = w2mean / wmean**2
+        # neff = ngal / (4 * np.pi * fsky) / deff
         bias = fsky * musq / dens
 
         # set metadata of array
@@ -453,11 +453,11 @@ async def __call__(
 
         # bias from measured variance, for E/B decomposition
         musq = var / wmean**2
-        dens = ngal / (4 * np.pi * fsky)  # should we include the factor of 2 here?
-        #variance = var / w2mean
-        #deff = w2mean / wmean**2
-        #neff = ngal / (2 * np.pi * fsky) / deff
-        bias = (1/2) * fsky * musq / dens 
+        dens = ngal / (4 * np.pi * fsky)
+        # variance = var / w2mean
+        # deff = w2mean / wmean**2
+        # neff = ngal / (2 * np.pi * fsky) / deff
+        bias = (1 / 2) * fsky * musq / dens
 
         # set metadata of array
         update_metadata(
@@ -559,10 +559,15 @@ async def __call__(
         wht /= wbar
 
         # bias from weights
-        bias = 4 * np.pi * fsky**2 * (w2mean / wmean**2) / ngal
+        musq = w2mean / wmean**2  # 1.0
+        dens = ngal / (4 * np.pi * fsky)
+        bias = fsky * musq / dens
+        # bias = 4 * np.pi * fsky**2 * (w2mean / wmean**2) / ngal
 
         # set metadata of array
-        update_metadata(wht, catalog, wbar=wbar, bias=bias)
+        update_metadata(
+            wht, catalog, wbar=wbar, musq=musq, dens=dens, fsky=fsky, bias=bias
+        )
 
         # return the weight map
         return wht

diff --git a/tests/test_fields.py b/tests/test_fields.py
@@ -306,7 +306,7 @@ def test_scalar_field(mapper, catalog):
     bias = (4 * np.pi / npix / npix) * v2
     wmean = v1 / (4.0 * npix)
     w2mean = v2 / (4.0 * npix)
-    variance = w2mean / wmean**2
+    musq = w2mean / wmean**2
 
     assert m.shape == (npix,)
     assert m.dtype.metadata == {
@@ -318,7 +318,7 @@ def test_scalar_field(mapper, catalog):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
-        "musq": pytest.approx(variance),
+        "musq": pytest.approx(musq),
         "dens": pytest.approx(npix / np.pi),
         "fsky": 1.0,
         "bias": pytest.approx(bias / wbar**2),
@@ -344,7 +344,7 @@ def test_complex_field(mapper, catalog):
     bias = (4 * np.pi / npix / npix) * v2 / 2
     wmean = v1 / (4.0 * npix)
     w2mean = v2 / (4.0 * npix)
-    variance = w2mean / wmean**2
+    musq = w2mean / wmean**2
 
     assert m.shape == (2, npix)
     assert m.dtype.metadata == {
@@ -358,7 +358,7 @@ def test_complex_field(mapper, catalog):
         "deconv": mapper.deconvolve,
         "fsky": 1.0,
         "dens": npix / np.pi,
-        "musq": pytest.approx(variance),
+        "musq": pytest.approx(musq),
         "bias": pytest.approx(bias / wbar**2),
     }
     np.testing.assert_array_almost_equal(m, 0)
@@ -377,6 +377,10 @@ def test_weights(mapper, catalog):
     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)
+    musq = w2mean / wmean**2
 
     assert m.shape == (12 * mapper.nside**2,)
     assert m.dtype.metadata == {
@@ -388,6 +392,9 @@ def test_weights(mapper, catalog):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
+        "musq": pytest.approx(musq),
+        "dens": pytest.approx(npix / np.pi),
+        "fsky": 1.0,
         "bias": pytest.approx(bias / wbar**2),
     }
     np.testing.assert_array_almost_equal(m, w / wbar)