From 685ab81c2ce752cd8fe77c8edf7237152ca4e912 Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Fri, 1 Mar 2024 06:25:12 -0600
Subject: [PATCH 1/7] ENH add fixed apodization for DES

---
 piff/psf.py | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 56 insertions(+), 1 deletion(-)

diff --git a/piff/psf.py b/piff/psf.py
index aa6f12ca..1bead0b9 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -20,10 +20,31 @@
 import fitsio
 import galsim
 import sys
+from numba import vectorize
 
 from .star import Star, StarData
 from .util import write_kwargs, read_kwargs
 
+
+@vectorize
+def _ap_kern_kern(x, m, h):
+    # cumulative triweight kernel
+    y = (x - m) / h + 3
+    if y < -3:
+        return 0
+    elif y > 3:
+        return 1
+    else:
+        val = (
+            -5 * y ** 7 / 69984
+            + 7 * y ** 5 / 2592
+            - 35 * y ** 3 / 864
+            + 35 * y / 96
+            + 1 / 2
+        )
+        return val
+
+
 class PSF(object):
     """The base class for describing a PSF model across a field of view.
 
@@ -99,7 +120,7 @@ def parseKwargs(cls, config_psf, logger=None):
         raise NotImplementedError("Derived classes must define the parseKwargs function")
 
     def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_size=48,
-             image=None, logger=None, **kwargs):
+             image=None, logger=None, apodize=True, **kwargs):
         r"""Draws an image of the PSF at a given location.
 
         The normal usage would be to specify (chipnum, x, y), in which case Piff will use the
@@ -161,6 +182,10 @@ def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_siz
                             [default: 48]
         :param image:       An existing image on which to draw, if desired. [default: None]
         :param logger:      A logger object for logging debug info. [default: None]
+        :param apodize:     Optional float giving the number of half light radii after
+                            which the profile is smoothy apodized to zero at the image edge.
+                            [default: True which uses a default transition width of 1 pixel and
+                            a default apodization radius of 4.25 pixels]
         :param \**kwargs:   Any additional properties required for the interpolation.
 
         :returns:           A GalSim Image of the PSF
@@ -201,6 +226,36 @@ def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_siz
 
         prof.drawImage(image, method=method, center=center)
 
+        if apodize:
+            xpix, ypix = image.get_pixel_centers()
+            dx = xpix - center[0]
+            dy = ypix - center[1]
+            r2 = dx**2 + dy**2
+
+            # this algorithm is adaptive but for DES we'll fix things
+            # to a constant value for stability
+            # _im = image.copy()
+            # _im.wcs = None
+            # _im.scale = 1.0
+            # hlr = _im.calculateHLR(center=galsim.PositionD(center))
+            # aprad = apodize * hlr
+            # msk_nonzero = _im.array != 0
+            # max_r = min(
+            #     np.abs(dx[(dx < 0) & msk_nonzero].min()),
+            #     np.abs(dx[(dx > 0) & msk_nonzero].max()),
+            #     np.abs(dy[(dy < 0) & msk_nonzero].min()),
+            #     np.abs(dy[(dy > 0) & msk_nonzero].max()),
+            # )
+            # apwidth = np.abs(hlr) / 2.355
+            # apwidth = min(max(apwidth, 0.5), 5.0)
+            # aprad = max(min(aprad, max_r - 6 * apwidth - 1), 2 * apwidth)
+
+            apwidth = 1.0
+            aprad = 4.25
+
+            apim = image._array * _ap_kern_kern(aprad, np.sqrt(r2), apwidth)
+            image._array = apim / np.sum(apim) * np.sum(image._array)
+
         return image
 
     def get_profile(self, x, y, chipnum=None, flux=1.0, logger=None, **kwargs):

From 54fba85bc50e3ebb0d9a9e5c607a3917095dbcac Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Fri, 1 Mar 2024 06:42:19 -0600
Subject: [PATCH 2/7] REF redo options a bit

---
 piff/psf.py | 52 +++++++++++++++++++++++++++-------------------------
 1 file changed, 27 insertions(+), 25 deletions(-)

diff --git a/piff/psf.py b/piff/psf.py
index 1bead0b9..5c3ff0c1 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -120,7 +120,7 @@ def parseKwargs(cls, config_psf, logger=None):
         raise NotImplementedError("Derived classes must define the parseKwargs function")
 
     def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_size=48,
-             image=None, logger=None, apodize=True, **kwargs):
+             image=None, logger=None, apodize=(1.0, 4.25), **kwargs):
         r"""Draws an image of the PSF at a given location.
 
         The normal usage would be to specify (chipnum, x, y), in which case Piff will use the
@@ -182,10 +182,11 @@ def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_siz
                             [default: 48]
         :param image:       An existing image on which to draw, if desired. [default: None]
         :param logger:      A logger object for logging debug info. [default: None]
-        :param apodize:     Optional float giving the number of half light radii after
-                            which the profile is smoothy apodized to zero at the image edge.
-                            [default: True which uses a default transition width of 1 pixel and
-                            a default apodization radius of 4.25 pixels]
+        :param apodize:     Optional parameter to set apodizatoon. If a float/int, gives the
+                            number of half light radii after which the profile is smoothy apodized
+                            to zero a width of ~2.55 half light radii. If a tuple/list, gives
+                            the apodization width and the apodization radius in pixels.
+                            [default: None, which means no apodization.]
         :param \**kwargs:   Any additional properties required for the interpolation.
 
         :returns:           A GalSim Image of the PSF
@@ -232,26 +233,27 @@ def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_siz
             dy = ypix - center[1]
             r2 = dx**2 + dy**2
 
-            # this algorithm is adaptive but for DES we'll fix things
-            # to a constant value for stability
-            # _im = image.copy()
-            # _im.wcs = None
-            # _im.scale = 1.0
-            # hlr = _im.calculateHLR(center=galsim.PositionD(center))
-            # aprad = apodize * hlr
-            # msk_nonzero = _im.array != 0
-            # max_r = min(
-            #     np.abs(dx[(dx < 0) & msk_nonzero].min()),
-            #     np.abs(dx[(dx > 0) & msk_nonzero].max()),
-            #     np.abs(dy[(dy < 0) & msk_nonzero].min()),
-            #     np.abs(dy[(dy > 0) & msk_nonzero].max()),
-            # )
-            # apwidth = np.abs(hlr) / 2.355
-            # apwidth = min(max(apwidth, 0.5), 5.0)
-            # aprad = max(min(aprad, max_r - 6 * apwidth - 1), 2 * apwidth)
-
-            apwidth = 1.0
-            aprad = 4.25
+            if isinstance(apodize, (tuple, list)):
+                apwidth, aprad = apodize
+            else:
+                wcs = image.wcs
+                try:
+                    image.wcs = None
+                    image.scale = 1.0
+                    hlr = image.calculateHLR(center=galsim.PositionD(center))
+                finally:
+                    image.wcs = wcs
+                aprad = apodize * hlr
+                msk_nonzero = image.array != 0
+                max_r = min(
+                    np.abs(dx[(dx < 0) & msk_nonzero].min()),
+                    np.abs(dx[(dx > 0) & msk_nonzero].max()),
+                    np.abs(dy[(dy < 0) & msk_nonzero].min()),
+                    np.abs(dy[(dy > 0) & msk_nonzero].max()),
+                )
+                apwidth = np.abs(hlr) / 2.355
+                apwidth = min(max(apwidth, 0.5), 5.0)
+                aprad = max(min(aprad, max_r - 6 * apwidth - 1), 2 * apwidth)
 
             apim = image._array * _ap_kern_kern(aprad, np.sqrt(r2), apwidth)
             image._array = apim / np.sum(apim) * np.sum(image._array)

From 139066ab8ca1ca23349bd0a3cc475c0c8182172a Mon Sep 17 00:00:00 2001
From: beckermr <becker.mr@gmail.com>
Date: Fri, 1 Mar 2024 06:45:46 -0600
Subject: [PATCH 3/7] REF no numba

---
 piff/psf.py | 27 ++++++++++++---------------
 1 file changed, 12 insertions(+), 15 deletions(-)

diff --git a/piff/psf.py b/piff/psf.py
index 5c3ff0c1..4a6496ad 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -20,29 +20,26 @@
 import fitsio
 import galsim
 import sys
-from numba import vectorize
 
 from .star import Star, StarData
 from .util import write_kwargs, read_kwargs
 
 
-@vectorize
 def _ap_kern_kern(x, m, h):
     # cumulative triweight kernel
     y = (x - m) / h + 3
-    if y < -3:
-        return 0
-    elif y > 3:
-        return 1
-    else:
-        val = (
-            -5 * y ** 7 / 69984
-            + 7 * y ** 5 / 2592
-            - 35 * y ** 3 / 864
-            + 35 * y / 96
-            + 1 / 2
-        )
-        return val
+    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 PSF(object):

From aba23c760f9e2f0d4cc87ef2941d8961cf438d25 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Sat, 2 Mar 2024 06:22:38 -0600
Subject: [PATCH 4/7] Update piff/psf.py

---
 piff/psf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/piff/psf.py b/piff/psf.py
index 4a6496ad..3850602a 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -117,7 +117,7 @@ def parseKwargs(cls, config_psf, logger=None):
         raise NotImplementedError("Derived classes must define the parseKwargs function")
 
     def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_size=48,
-             image=None, logger=None, apodize=(1.0, 4.25), **kwargs):
+             image=None, logger=None, apodize=None, **kwargs):
         r"""Draws an image of the PSF at a given location.
 
         The normal usage would be to specify (chipnum, x, y), in which case Piff will use the

From 8e645842068ee49c3d00facd6c869436cd56eb96 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:48:17 -0500
Subject: [PATCH 5/7] Update piff/psf.py

---
 piff/psf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/piff/psf.py b/piff/psf.py
index 3850602a..4a6496ad 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -117,7 +117,7 @@ def parseKwargs(cls, config_psf, logger=None):
         raise NotImplementedError("Derived classes must define the parseKwargs function")
 
     def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_size=48,
-             image=None, logger=None, apodize=None, **kwargs):
+             image=None, logger=None, apodize=(1.0, 4.25), **kwargs):
         r"""Draws an image of the PSF at a given location.
 
         The normal usage would be to specify (chipnum, x, y), in which case Piff will use the

From c76660a2db16fac4b998b3f8d45fde0e12f9e4c2 Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:48:59 -0500
Subject: [PATCH 6/7] Update piff/psf.py

---
 piff/psf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/piff/psf.py b/piff/psf.py
index 4a6496ad..a6e935cf 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -183,7 +183,7 @@ def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_siz
                             number of half light radii after which the profile is smoothy apodized
                             to zero a width of ~2.55 half light radii. If a tuple/list, gives
                             the apodization width and the apodization radius in pixels.
-                            [default: None, which means no apodization.]
+                            [default: (1.0, 4.25), which means a width of 1 pixel and radius of 4.25 pixels.]
         :param \**kwargs:   Any additional properties required for the interpolation.
 
         :returns:           A GalSim Image of the PSF

From cec733c7979d5854571cb91b3b7fde7807f0494e Mon Sep 17 00:00:00 2001
From: "Matthew R. Becker" <beckermr@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:51:31 -0500
Subject: [PATCH 7/7] Update _version.py

---
 piff/_version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/piff/_version.py b/piff/_version.py
index 841c3f9e..c1ba101f 100644
--- a/piff/_version.py
+++ b/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'
+__version__ = '1.3.3.1'
 __version_info__ = tuple(map(int, __version__.split('.')))