ENH apodize for DES

.github/workflows/ci.yml

@@ -16,6 +16,7 @@ jobs:
         runs-on: ${{ matrix.os }}
 
         strategy:
+            fail-fast: false
             matrix:
                 # First all python versions in basic linux
                 os: [ ubuntu-latest ]
@@ -81,13 +82,13 @@ jobs:
                 python -m pip install -U pip
 
                 # Do these first to clarify potential conflicts
-                pip install -U setuptools numpy
+                pip install -U setuptools "numpy<2"
 
                 # Standard dependencies
                 pip install -U -r requirements.txt
 
                 # Extra packages needed for testing
-                pip install -U nose coverage pytest nbval ipykernel
+                pip install -U nose coverage "pytest<8" nbval ipykernel
 
             - name: Install Pixmappy (not on pip)
               run: |

piff/_version.py

@@ -12,5 +12,5 @@
 #    this list of conditions and the disclaimer given in the documentation
 #    and/or other materials provided with the distribution.
 
-__version__ = '1.3.3'
+__version__ = '1.3.3.1'
 __version_info__ = tuple(map(int, __version__.split('.')))

piff/pixelgrid.py

@@ -25,6 +25,31 @@
 from .model import Model
 from .star import Star, StarData, StarFit
 
+APODIZE_PARAMS = (1.0 * 0.263, 4.25 * 0.263)
+
+
+def set_apodize_params(pars):
+    global APODIZE_PARAMS
+    APODIZE_PARAMS = pars
+
+
+def _ap_kern_kern(x, m, h):
+    # cumulative triweight kernel
+    y = (x - m) / h + 3
+    apval = np.zeros_like(m)
+    msk = y > 3
+    apval[msk] = 1
+    msk = (y > -3) & (~msk)
+    apval[msk] = (
+        -5 * y[msk] ** 7 / 69984
+        + 7 * y[msk] ** 5 / 2592
+        - 35 * y[msk] ** 3 / 864
+        + 35 * y[msk] / 96
+        + 1 / 2
+    )
+    return apval
+
+
 class PixelGrid(Model):
     """A PSF modeled as interpolation between a grid of points.
 
@@ -445,6 +470,21 @@ def getProfile(self, params):
         :returns: a galsim.GSObject instance
         """
         im = galsim.Image(params.reshape(self.size,self.size), scale=self.scale)
+
+        if APODIZE_PARAMS is not None:
+            xpix, ypix = im.get_pixel_centers()
+            # use_true_center = False below
+            dx = xpix - im.center.x
+            dy = ypix - im.center.y
+            r2 = dx**2 + dy**2
+
+            apwidth, aprad = APODIZE_PARAMS  # in arcsec
+            _apwidth = apwidth / self.scale  # convert to pixels
+            _aprad = aprad / self.scale  # convert to pixels
+
+            apim = im._array * _ap_kern_kern(_aprad, np.sqrt(r2), _apwidth)
+            im._array = apim / np.sum(apim) * np.sum(im._array)
+
         return galsim.InterpolatedImage(im, x_interpolant=self.interp,
                                         use_true_center=False, flux=1.)
 

requirements.txt

@@ -1,10 +1,10 @@
-numpy>=1.17
+numpy>=1.17,<2
 scipy>=1.2
 pyyaml>=5.1
 lsstdesc.coord>=1.2
 treecorr>=4.3.1
 fitsio>=1.0
 matplotlib>=3.3
-galsim>=2.3
+galsim>=2.3,<2.5
 treegp>=0.6
 threadpoolctl>=3.1

tests/conftest.py

@@ -0,0 +1,4 @@
+# turn off apodization
+import piff.pixelgrid
+
+piff.pixelgrid.set_apodize_params(None)

tests/test_wcs.py

@@ -747,6 +747,53 @@ def test_des_wcs():
     np.testing.assert_allclose(wcs3.toWorld(im.center).y, wcs1.toWorld(im.center).y,
                                rtol=0.04)
 
+@timer
+def test_newdes_apodize():
+    # This is a DES Y6 PSF file made by Robert Gruendl using python 2, so
+    # check that this also works correctly.
+    try:
+        import pixmappy
+    except ImportError:
+        print('pixmappy not installed.  Skipping test_newdes()')
+        return
+    # Also make sure pixmappy is recent enough to work.
+    if 'exposure_file' not in pixmappy.GalSimWCS._opt_params:
+        print('pixmappy not recent enough version.  Skipping test_newdes()')
+        return
+
+    import piff
+    import piff.pixelgrid
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_newdes.log')
+
+    fname = os.path.join('input', 'D00232418_i_c19_r5006p01_piff-model.fits')
+    with warnings.catch_warnings():
+        # This file was written with GalSim 2.1, and now raises a deprecation warning for 2.2.
+        warnings.simplefilter("ignore", galsim.GalSimDeprecationWarning)
+        warnings.simplefilter("ignore", DeprecationWarning)
+        psf = piff.PSF.read(fname, logger=logger)
+
+    ims = []
+    for appars in [None, (1.0 * 0.263, 4.25 * 0.263)]:
+        piff.pixelgrid.set_apodize_params(appars)
+        ims.append(psf.draw(x=103.3, y=592.0, logger=logger))
+
+    print('sum = ',ims[1].array.sum())
+    assert not np.allclose(ims[0].array, ims[1].array)
+    assert np.allclose(ims[1].array[0, :], 0, rtol=1.e-2)
+    assert np.allclose(ims[1].array[-1, :], 0, rtol=1.e-2)
+    assert np.allclose(ims[1].array[:, 0], 0, rtol=1.e-2)
+    assert np.allclose(ims[1].array[:, -1], 0, rtol=1.e-2)
+    assert ims[1].array.sum() > 0
+    np.testing.assert_allclose(
+        ims[0].array[23:26,22:25] / ims[0].array[23:26,22:25].sum(),
+        ims[1].array[23:26,22:25] / ims[1].array[23:26,22:25].sum(),
+        rtol=1.e-5,
+    )
+
 if __name__ == '__main__':
     #import cProfile, pstats
     #pr = cProfile.Profile()