Add psf_fwhm (median seeing) for make_ps_psf. Add parameters for make_rand_psf

descwl_shear_sims/psfs/ps_psf.py

@@ -1,7 +1,7 @@
 """
 Copied from https://github.com/beckermr/metadetect-sims under BSD
 """
-from ..constants import SCALE
+from ..constants import SCALE, FIXED_PSF_FWHM
 
 import numpy as np
 import galsim
@@ -10,7 +10,14 @@
 import galsim.utilities
 
 
-def make_ps_psf(*, rng, dim, pixel_scale=SCALE, variation_factor=1):
+def make_ps_psf(
+    *,
+    rng,
+    dim,
+    pixel_scale=SCALE,
+    psf_fwhm=FIXED_PSF_FWHM,
+    variation_factor=1,
+):
     """
     get a power spectrum psf
 
@@ -22,6 +29,8 @@ def make_ps_psf(*, rng, dim, pixel_scale=SCALE, variation_factor=1):
         Dimensions of image
     pixel_scale: float
         pixel scale
+    psf_fwhm: float
+       Median FWHM of PSF in units of arcsec
     variation_factor : float, optional
         This factor is used internally to scale the overall variance in the
         PSF shape power spectra and the change in the PSF size across the
@@ -37,6 +46,7 @@ def make_ps_psf(*, rng, dim, pixel_scale=SCALE, variation_factor=1):
         im_width=dim,
         buff=dim/2,
         scale=pixel_scale,
+        median_seeing=psf_fwhm,
         variation_factor=variation_factor,
     )
 

descwl_shear_sims/psfs/rand_psf.py

@@ -11,7 +11,16 @@
 )
 
 
-def make_rand_psf(psf_type, rng):
+def make_rand_psf(
+    psf_type,
+    rng,
+    psf_fwhm_mean=RAND_PSF_FWHM_MEAN,
+    psf_fwhm_std=RAND_PSF_FWHM_STD,
+    psf_fwhm_min=RAND_PSF_FWHM_MIN,
+    psf_fwhm_max=RAND_PSF_FWHM_MAX,
+    psf_e_std=RAND_PSF_E_STD,
+    psf_e_max=RAND_PSF_E_MAX,
+):
     """
     A simple PSF object with random FWHM and shape
 
@@ -25,8 +34,18 @@ def make_rand_psf(psf_type, rng):
     Gaussian or Moffat
     """
 
-    fwhm = _get_fwhm(rng)
-    e1, e2 = _get_e1e2(rng)
+    fwhm = _get_fwhm(
+        rng,
+        fwhm_mean=psf_fwhm_mean,
+        fwhm_std=psf_fwhm_std,
+        fwhm_min=psf_fwhm_min,
+        fwhm_max=psf_fwhm_max,
+    )
+    e1, e2 = _get_e1e2(
+        rng,
+        e_std=psf_e_std,
+        e_max=psf_e_max,
+    )
 
     if psf_type == "gauss":
         psf = galsim.Gaussian(fwhm=fwhm)
@@ -40,28 +59,28 @@ def make_rand_psf(psf_type, rng):
     return psf
 
 
-def _get_fwhm(rng):
+def _get_fwhm(rng, fwhm_mean, fwhm_std, fwhm_min, fwhm_max):
     ln_mean = log(
-        RAND_PSF_FWHM_MEAN**2 / sqrt(RAND_PSF_FWHM_MEAN**2 + RAND_PSF_FWHM_STD**2)
+        fwhm_mean**2 / sqrt(fwhm_mean**2 + fwhm_std**2)
     )  # noqa
-    ln_sigma = sqrt(log(1+(RAND_PSF_FWHM_STD/RAND_PSF_FWHM_MEAN)**2))
+    ln_sigma = sqrt(log(1+(fwhm_std/fwhm_mean)**2))
 
     while True:
         fwhm = rng.lognormal(
             mean=ln_mean,
             sigma=ln_sigma,
         )
-        if RAND_PSF_FWHM_MIN < fwhm < RAND_PSF_FWHM_MAX:
+        if fwhm_min < fwhm < fwhm_max:
             break
 
     return fwhm
 
 
-def _get_e1e2(rng):
+def _get_e1e2(rng, e_std, e_max):
     while True:
-        e1, e2 = rng.normal(scale=RAND_PSF_E_STD, size=2)
+        e1, e2 = rng.normal(scale=e_std, size=2)
         e = sqrt(e1**2 + e2**2)
-        if e < RAND_PSF_E_MAX:
+        if e < e_max:
             break
 
     return e1, e2

descwl_shear_sims/sim.py

@@ -27,6 +27,10 @@
     "gal_type": "fixed",  # note "exp" also means "fixed" for back compat
     "psf_type": "gauss",
     "psf_dim": 51,
+    # constant FWHM value for fixed PSF
+    # mean for random PSF
+    # median for PS psf
+    "psf_fwhm": 0.8,
     "psf_variation_factor": 1,  # for power spectrum psf
     # randomize the psf fwhm and shape for each trial.  PSF is still same
     # for all epochs/bands

descwl_shear_sims/tests/test_config.py

@@ -8,9 +8,10 @@ def test_config_smoke():
 
 
 def test_config_input():
-    cin = {'psf_type': 'moffat'}
+    cin = {'psf_type': 'moffat', 'psf_fwhm': 0.9}
     config = get_sim_config(config=cin)
     assert config['psf_type'] == cin['psf_type']
+    assert config['psf_fwhm'] == cin['psf_fwhm']
 
 
 def test_config_pure():

descwl_shear_sims/tests/test_ps_psf.py

@@ -34,7 +34,17 @@ def test_ps_psf_smoke():
     psf_im = psf.drawImage(scale=PIXEL_SCALE)
     assert psf_im.calculateFWHM() > 0.5
 
-    ps = make_ps_psf(rng=rng, dim=120)
+    fwhm = 0.9
+    varfac = 1.2
+    ps = make_ps_psf(
+        rng=rng,
+        dim=120,
+        psf_fwhm=fwhm,
+        variation_factor=varfac,
+    )
+    assert ps._median_seeing == 0.9
+    assert ps._variation_factor == varfac
+
     psf = ps.getPSF(galsim.PositionD(x=10, y=10))
     assert isinstance(psf, galsim.GSObject)
     psf_im = psf.drawImage(scale=PIXEL_SCALE)

descwl_shear_sims/tests/test_rand_psf.py

@@ -4,14 +4,15 @@
 
 import numpy as np
 import galsim
+import pytest
 
 # import pytest
 import lsst.geom as geom
 import lsst.afw.image as afw_image
 
 from ..psfs import make_rand_psf, make_dm_psf
 from ._wcs import make_sim_wcs, SCALE
-from ..constants import RAND_PSF_FWHM_STD, RAND_PSF_E_STD
+from ..constants import RAND_PSF_FWHM_MEAN, RAND_PSF_FWHM_STD, RAND_PSF_E_STD
 
 
 def _get_fwhm_g1g2(psf_im):
@@ -27,18 +28,34 @@ def test_rand_psf_smoke():
     rng = np.random.RandomState(seed=10)
     psf = make_rand_psf(
         psf_type='gauss',
+        psf_fwhm_mean=0.85,
+        psf_fwhm_std=0.12,
+        psf_fwhm_min=0.6,
+        psf_fwhm_max=1.5,
+        psf_e_std=0.03,
+        psf_e_max=0.2,
         rng=rng,
     )
     assert isinstance(psf, galsim.GSObject)
 
 
-def test_rand_dmpsf_smoke():
+@pytest.mark.parametrize('defaults', [True, False])
+def test_rand_dmpsf_smoke(defaults):
+
+    if defaults:
+        fwhm_mean = RAND_PSF_FWHM_MEAN
+        fwhm_std = RAND_PSF_FWHM_STD
+        e_std = RAND_PSF_E_STD
+    else:
+        fwhm_mean = 0.85
+        fwhm_std = 0.12
+        e_std = 0.03
 
     dim = 20
     psf_dim = 31
     rng = np.random.RandomState(seed=999)
 
-    ntrial = 100
+    ntrial = 500
     fwhmvals = np.zeros(ntrial)
     e1vals = np.zeros(ntrial)
     e2vals = np.zeros(ntrial)
@@ -50,6 +67,13 @@ def test_rand_dmpsf_smoke():
             psf_type='gauss',
             rng=rng,
         )
+        gspsf = make_rand_psf(
+            psf_type='gauss',
+            psf_fwhm_mean=fwhm_mean,
+            psf_fwhm_std=fwhm_std,
+            psf_e_std=e_std,
+            rng=rng,
+        )
 
         galsim_wcs = make_sim_wcs(dim)
         # dm_wcs = make_dm_wcs(galsim_wcs)
@@ -74,9 +98,18 @@ def test_rand_dmpsf_smoke():
         e1vals[i] = e1
         e2vals[i] = e2
 
-    fwhmstd = fwhmvals.std()
-    e1std = e1vals.std()
-    e2std = e2vals.std()
-    assert fwhmstd/RAND_PSF_FWHM_STD - 1 < 0.05
-    assert e1std/RAND_PSF_E_STD - 1 < 0.05
-    assert e2std/RAND_PSF_E_STD - 1 < 0.05
+    afwhm_mean = fwhmvals.mean()
+    afwhm_std = fwhmvals.std()
+
+    ae1std = e1vals.std()
+    ae2std = e2vals.std()
+
+    fwhm_fracdiff = afwhm_mean / fwhm_mean - 1
+    fwhm_std_fracdiff = afwhm_std / fwhm_std - 1
+    e1std_fracdiff = ae1std / e_std - 1
+    e2std_fracdiff = ae2std / e_std - 1
+
+    assert fwhm_fracdiff < 0.05
+    assert fwhm_std_fracdiff < 0.05
+    assert e1std_fracdiff < 0.05
+    assert e2std_fracdiff < 0.05