From 573d557a7d85408d3b19456f7ab844f995a693ab Mon Sep 17 00:00:00 2001
From: jaimerzp <jaimerz011235813@gmail.com>
Date: Wed, 27 Nov 2024 15:20:45 +0000
Subject: [PATCH] working weghts

---
 heracles/fields.py   | 27 ++++++++++++++++-----------
 tests/test_fields.py | 15 +++++++++++----
 2 files changed, 27 insertions(+), 15 deletions(-)

diff --git a/heracles/fields.py b/heracles/fields.py
index e9397eb..872ef97 100644
--- 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
index 6d96a7e..db5a4e6 100644
--- 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)