remove magnification and rerun experiments (#79)

* remove magnification * calculate std of ellipticity * move notebooks around * rerun exp30 * rerun experiments no magnification
LSSTDESC · Jan 14, 2025 · 861b34f · 861b34f
1 parent 070f885
commit 861b34f
Show file tree

Hide file tree

Showing 23 changed files with 92 additions and 42 deletions.
diff --git a/bpd/draw.py b/bpd/draw.py
@@ -44,11 +44,7 @@ def draw_gaussian_galsim(
 ):
     gal = galsim.Gaussian(flux=f, half_light_radius=hlr)
     gal = gal.shear(g1=e1, g2=e2)  # intrinsic ellipticity
-
-    # the correct weak lensing effect includes magnification even if kappa=0!
-    # see: https://galsim-developers.github.io/GalSim/_build/html/shear.html
-    mu = (1 - g1**2 - g2**2) ** -1  # convergence kappa = 0
-    gal = gal.lens(g1=g1, g2=g2, mu=mu)
+    gal = gal.shear(g1=g1, g2=g2)
 
     psf = galsim.Gaussian(flux=1.0, half_light_radius=psf_hlr)
     gal_conv = galsim.Convolve([gal, psf])

diff --git a/experiments/exp30/figs/contours.pdf b/experiments/exp30/figs/contours.pdf
diff --git a/experiments/exp30/figs/hists.pdf b/experiments/exp30/figs/hists.pdf
diff --git a/experiments/exp30/figs/scatter_shapes.pdf b/experiments/exp30/figs/scatter_shapes.pdf
diff --git a/experiments/exp30/figs/traces.pdf b/experiments/exp30/figs/traces.pdf
diff --git a/experiments/exp31/figs/contours.pdf b/experiments/exp31/figs/contours.pdf
diff --git a/experiments/exp31/figs/hists.pdf b/experiments/exp31/figs/hists.pdf
diff --git a/experiments/exp31/figs/scatter_dxdy.pdf b/experiments/exp31/figs/scatter_dxdy.pdf
diff --git a/experiments/exp31/figs/scatter_shapes.pdf b/experiments/exp31/figs/scatter_shapes.pdf
diff --git a/experiments/exp31/figs/traces.pdf b/experiments/exp31/figs/traces.pdf
diff --git a/experiments/exp32/figs/contours.pdf b/experiments/exp32/figs/contours.pdf
diff --git a/experiments/exp32/figs/hists.pdf b/experiments/exp32/figs/hists.pdf
diff --git a/experiments/exp32/figs/hists_flux_hlr.pdf b/experiments/exp32/figs/hists_flux_hlr.pdf
diff --git a/experiments/exp32/figs/scatter_dxdy.pdf b/experiments/exp32/figs/scatter_dxdy.pdf
diff --git a/experiments/exp32/figs/scatter_shapes.pdf b/experiments/exp32/figs/scatter_shapes.pdf
diff --git a/experiments/exp32/figs/traces.pdf b/experiments/exp32/figs/traces.pdf
diff --git a/experiments/exp32/get_shear.py b/experiments/exp32/get_shear.py
@@ -42,16 +42,9 @@ def _logprior(
 
     prior = jnp.array(0.0)
 
-    # both flux and HLR get magnified by WL so we need to transform them
-    mu = (1 - g[0] ** 2 - g[1] ** 2) ** -1
-    lf_unsheared = lf - jnp.log(mu)
-    lhlr_unsheared = lhlr - 0.5 * jnp.log(mu)
-
     # we also pickup a jacobian in the corresponding probability densities
-    prior += stats.norm.logpdf(lf_unsheared, loc=mean_logflux, scale=sigma_logflux)
-    prior -= jnp.log(mu)  # jacobian (flux)
-    prior += stats.norm.logpdf(lhlr_unsheared, loc=mean_loghlr, scale=sigma_loghlr)
-    prior -= jnp.log(mu) * 0.5  # jacobian (hlr)
+    prior += stats.norm.logpdf(lf, loc=mean_logflux, scale=sigma_logflux)
+    prior += stats.norm.logpdf(lhlr, loc=mean_loghlr, scale=sigma_loghlr)
 
     # elliptcity
     prior += true_ellip_logprior(e_post, g, sigma_e=sigma_e)

diff --git a/notebooks/bulge_disk_exp1.ipynb → notebooks/future/bulge_disk_exp1.ipynb b/notebooks/bulge_disk_exp1.ipynb → notebooks/future/bulge_disk_exp1.ipynb
diff --git a/notebooks/spergel-speed1.ipynb → notebooks/future/spergel-speed1.ipynb b/notebooks/spergel-speed1.ipynb → notebooks/future/spergel-speed1.ipynb
diff --git a/notebooks/test-shear-inf-with-prior1.ipynb → ...s/future/test-shear-inf-with-prior1.ipynb b/notebooks/test-shear-inf-with-prior1.ipynb → ...s/future/test-shear-inf-with-prior1.ipynb
diff --git a/notebooks/old/back-envelop-ellip-mag1.ipynb b/notebooks/old/back-envelop-ellip-mag1.ipynb
@@ -0,0 +1,83 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "cca60799-f4c2-4395-ab9d-5461a2419123",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "from bpd.prior import sample_noisy_ellipticities_unclipped\n",
+    "from jax import random\n",
+    "import jax.numpy as jnp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "38140b10-b976-46e8-9d1e-da98b867acdf",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "rng_key = random.key(42)\n",
+    "_, e, _ = sample_noisy_ellipticities_unclipped(rng_key, g=jnp.array([0.02, 0.0]), sigma_m=1e-4, sigma_e=1e-3, n=10000)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "3fcea711-9436-4509-a6cb-52e61ae6a020",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Array(3.13401326e-05, dtype=float64)"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "e[:, 0].std() / jnp.sqrt(1000)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f0b78fb1-9fe8-4723-b31f-bf0cb5bc95e3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "bpd_gpu3",
+   "language": "python",
+   "name": "bpd_gpu3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/notebooks/magnitfy-test1.ipynb → notebooks/old/magnitfy-test1.ipynb b/notebooks/magnitfy-test1.ipynb → notebooks/old/magnitfy-test1.ipynb
diff --git a/tests/test_shear_trans.py b/tests/test_shear_trans.py
@@ -28,7 +28,6 @@ def _draw_gaussian(
     y: float,
     g1: float,
     g2: float,
-    mu: float,
     slen: int,
     fft_size: int,
     psf_hlr: float = 0.7,
@@ -38,7 +37,7 @@ def _draw_gaussian(
 
     gal = xgalsim.Gaussian(flux=f, half_light_radius=hlr)
     gal = gal.shear(g1=e1, g2=e2)
-    gal = gal.lens(g1=g1, g2=g2, mu=mu)
+    gal = gal.shear(g1=g1, g2=g2)
 
     psf = xgalsim.Gaussian(flux=1.0, half_light_radius=psf_hlr)
     gal_conv = xgalsim.Convolve([gal, psf]).withGSParams(gsparams)
@@ -96,38 +95,17 @@ def test_image_shear_commute():
                     g = jnp.array([g1, g2])
                     (e1_p, e2_p) = scalar_shear_transformation(e, g)
 
-                    # magnification
-                    mu = (1 - g1**2 - g2**2) ** -1
-                    f_p = f * mu
-                    hlr_p = hlr * mu ** (0.5)
-
-                    im1 = draw_jitted(
-                        f=f, hlr=hlr, e1=e1, e2=e2, g1=0, g2=0, mu=1.0, x=x, y=y
-                    )
+                    im1 = draw_jitted(f=f, hlr=hlr, e1=e1, e2=e2, g1=0, g2=0, x=x, y=y)
                     im2 = draw_jitted(
-                        f=f, hlr=hlr, e1=e1, e2=e2, g1=g1, g2=g2, mu=1.0, x=x, y=y
+                        f=f, hlr=hlr, e1=e1, e2=e2, g1=g1, g2=g2, x=x, y=y
                     )
                     im3 = draw_jitted(
-                        f=f, hlr=hlr, e1=e1_p, e2=e2_p, g1=0, g2=0, mu=1.0, x=x, y=y
-                    )
-                    im4 = draw_jitted(
-                        f=f, hlr=hlr, e1=e1, e2=e2, g1=g1, g2=g2, mu=mu, x=x, y=y
-                    )
-                    im5 = draw_jitted(
-                        f=f_p, hlr=hlr_p, e1=e1_p, e2=e2_p, g1=0, g2=0, mu=1.0, x=x, y=y
+                        f=f, hlr=hlr, e1=e1_p, e2=e2_p, g1=0, g2=0, x=x, y=y
                     )
 
                     # rtol is 0 because image contains lots of 0s
                     np.testing.assert_allclose(im2, im3, rtol=0, atol=1e-10)
-                    np.testing.assert_allclose(im4, im5, rtol=0, atol=5e-7)
 
                     if not (g1 == 0 and g2 == 0):
-                        assert not np.allclose(im1, im2, rtol=0, atol=1e-6)
-                        assert not np.allclose(im1, im3, rtol=0, atol=1e-6)
-                        assert not np.allclose(im1, im4, rtol=0, atol=1e-6)
-                        assert not np.allclose(im1, im5, rtol=0, atol=1e-6)
-
-                        assert not np.allclose(im2, im4, rtol=0, atol=1e-6)
-                        assert not np.allclose(im2, im5, rtol=0, atol=1e-6)
-                        assert not np.allclose(im3, im5, rtol=0, atol=1e-6)
-                        assert not np.allclose(im3, im5, rtol=0, atol=1e-6)
+                        assert not np.allclose(im1, im2, rtol=0, atol=1e-10)
+                        assert not np.allclose(im1, im3, rtol=0, atol=1e-10)