gh-212: metadata for bias computation

heracles/fields.py

@@ -305,10 +305,14 @@ async def __call__(
             del vis
 
         # compute bias of positions, including weight variance
-        bias = ngal / (4 * np.pi) * mapper.area**2 * (w2mean / nbar**2)
+        variance = w2mean / wmean**2
+        neff = (npix * nbar) / (4 * np.pi * wmean)
+        bias = fsky * variance / neff
 
         # set metadata of array
-        update_metadata(pos, catalog, nbar=nbar, bias=bias)
+        update_metadata(
+            pos, catalog, nbar=nbar, variance=variance, neff=neff, fsky=fsky, bias=bias
+        )
 
         # return the position map
         return pos
@@ -371,10 +375,14 @@ async def __call__(
         val /= wbar
 
         # compute bias from variance (per object)
-        bias = 4 * np.pi * fsky**2 * (var / wmean**2) / ngal
+        variance = var / wmean**2
+        neff = ngal / (4 * np.pi * fsky)
+        bias = fsky * variance / neff
 
         # set metadata of array
-        update_metadata(val, catalog, wbar=wbar, bias=bias)
+        update_metadata(
+            val, catalog, wbar=wbar, variance=variance, neff=neff, fsky=fsky, bias=bias
+        )
 
         # return the value map
         return val
@@ -441,10 +449,14 @@ async def __call__(
         val /= wbar
 
         # bias from measured variance, for E/B decomposition
-        bias = 2 * np.pi * fsky**2 * (var / wmean**2) / ngal
+        variance = var / wmean**2
+        neff = ngal / (2 * np.pi * fsky)  # should we include the factor of 2 here?
+        bias = fsky * variance / neff
 
         # set metadata of array
-        update_metadata(val, catalog, wbar=wbar, bias=bias)
+        update_metadata(
+            val, catalog, wbar=wbar, variance=variance, neff=neff, fsky=fsky, bias=bias
+        )
 
         # return the shear map
         return val

tests/test_fields.py

@@ -204,6 +204,9 @@ def test_positions(mapper, catalog, vmap):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
+        "variance": 1.0,
+        "neff": npix / np.pi,
+        "fsky": 1.0,
         "bias": pytest.approx(bias / nbar**2),
     }
     np.testing.assert_array_equal(m, 0)
@@ -223,6 +226,9 @@ def test_positions(mapper, catalog, vmap):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
+        "variance": 1.0,
+        "neff": npix / np.pi,
+        "fsky": 1.0,
         "bias": pytest.approx(bias / nbar**2),
     }
     np.testing.assert_array_equal(m, 1.0)
@@ -246,6 +252,9 @@ def test_positions(mapper, catalog, vmap):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
+        "variance": 1.0,
+        "neff": npix / (np.pi * catalog.fsky),
+        "fsky": catalog.fsky,
         "bias": pytest.approx(bias / nbar**2),
     }
 
@@ -264,6 +273,9 @@ def test_positions(mapper, catalog, vmap):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
+        "variance": 1.0,
+        "neff": npix / (np.pi * catalog.fsky),
+        "fsky": catalog.fsky,
         "bias": pytest.approx(bias / nbar**2),
     }
 
@@ -274,7 +286,7 @@ 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)
+    assert m.dtype.metadata["bias"] == pytest.approx(2 * bias / (2 * nbar) ** 2)
 
 
 def test_scalar_field(mapper, catalog):
@@ -290,7 +302,11 @@ def test_scalar_field(mapper, catalog):
     v2 = ((w * v) ** 2).sum()
     w = w.reshape(w.size // 4, 4).sum(axis=-1)
     wbar = w.mean()
+    v1 = w.sum()
     bias = (4 * np.pi / npix / npix) * v2
+    wmean = v1 / (4.0 * npix)
+    w2mean = v2 / (4.0 * npix)
+    variance = w2mean / wmean**2
 
     assert m.shape == (npix,)
     assert m.dtype.metadata == {
@@ -302,6 +318,9 @@ def test_scalar_field(mapper, catalog):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
+        "variance": pytest.approx(variance),
+        "neff": npix / np.pi,
+        "fsky": 1.0,
         "bias": pytest.approx(bias / wbar**2),
     }
     np.testing.assert_array_almost_equal(m, 0)
@@ -321,7 +340,11 @@ def test_complex_field(mapper, catalog):
     v2 = ((w * re) ** 2 + (w * im) ** 2).sum()
     w = w.reshape(w.size // 4, 4).sum(axis=-1)
     wbar = w.mean()
+    v1 = w.sum()
     bias = (4 * np.pi / npix / npix) * v2 / 2
+    wmean = v1 / (4.0 * npix)
+    w2mean = v2 / (4.0 * npix)
+    variance = w2mean / wmean**2
 
     assert m.shape == (2, npix)
     assert m.dtype.metadata == {
@@ -333,6 +356,9 @@ def test_complex_field(mapper, catalog):
         "nside": mapper.nside,
         "lmax": mapper.lmax,
         "deconv": mapper.deconvolve,
+        "fsky": 1.0,
+        "neff": 2 * npix / np.pi,
+        "variance": pytest.approx(variance),
         "bias": pytest.approx(bias / wbar**2),
     }
     np.testing.assert_array_almost_equal(m, 0)