From a05907bc792cf01709f2b453b5e566f73e1a99c4 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Tue, 1 Oct 2024 17:04:51 +0530
Subject: [PATCH 01/14] Created new PSFExPSF class

---
 photutils/psf/image_models.py | 78 +++++++++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 273780475..d37616c2b 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -331,6 +331,84 @@ def evaluate(self, x, y, flux, x_0, y_0):
 
         return evaluated_model
 
+class PSFExPSF(ImagePSF):
+    """
+    A fittable model for PSF output from PSFEx.
+
+    The PSFEx PSF is a grid of 2D models where each 2D element represents
+    the contribution of a polynomial. The polymial can be constructed from
+    any FITS header parameter, but the most straightforward example would
+    be the x and y positions of the CCD image.
+
+    The model has three model parameters: an image intensity scaling
+    factor (``flux``) which is applied to the input image, and two
+    positional parameters (``x_0`` and ``y_0``) indicating the location
+    of a feature in the coordinate grid on which the model is evaluated.
+
+    Parameters
+    ----------
+    data : 3D `~numpy.ndarray`
+        Array containing the grid of reference PSF arrays. The length of 
+        the x and y axes must both be at least 4. All elements of the 
+        input image data must be finite. By default, the PSF peak is 
+        assumed to be located at the center of the input image (see the ``origin``
+        keyword).
+
+    flux : float, optional
+        The total flux of the source, assuming the input image
+        was properly normalized.
+
+    x_0, y_0 : float
+        The x and y positions of a feature in the image in the output
+        coordinate grid on which the model is evaluated. Typically, this
+        refers to the position of the PSF peak, which is assumed to be
+        located at the center of the input image (see the ``origin``
+        keyword).
+
+    origin : tuple of 2 float or None, optional
+        The ``(x, y)`` coordinate with respect to the input image data
+        array that represents the reference pixel of the input data.
+
+        The reference ``origin`` pixel will be placed at the model
+        ``x_0`` and ``y_0`` coordinates in the output coordinate system
+        on which the model is evaluated.
+
+        Most typically, the input PSF should be centered in the input
+        image, and thus the origin should be set to the central pixel of
+        the ``data`` array.
+
+        If the origin is set to `None`, then the origin will be set to
+        the center of the ``data`` array (``(npix - 1) / 2.0``).
+
+    oversampling : int or array_like (int), optional
+        The integer oversampling factor(s) of the input PSF image. If
+        ``oversampling`` is a scalar then it will be used for both axes.
+        If ``oversampling`` has two elements, they must be in ``(y, x)``
+        order.
+
+    fill_value : float, optional
+        The value to use for points outside of the input pixel grid.
+        The default is 0.0.
+
+    **kwargs : dict, optional
+        Additional optional keyword arguments to be passed to the
+        `astropy.modeling.Model` base class.
+
+    See Also
+    --------
+    GriddedPSFModel : A model for a grid of ePSF models.    
+
+    """
+    flux = Parameter(default=1,
+                     description='Intensity scaling factor of the image.')
+    x_0 = Parameter(default=0,
+                    description=('Position of a feature in the image along '
+                                 'the x axis'))
+    y_0 = Parameter(default=0,
+                    description=('Position of a feature in the image along '
+                                 'the y axis'))
+
+    
 
 @deprecated('2.0.0', alternative='`ImagePSF`')
 class FittableImageModel(Fittable2DModel):

From da99af13334ca149b5dbb2b16e37b847f437e49a Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Thu, 3 Oct 2024 17:03:02 +0530
Subject: [PATCH 02/14] Update 1 to PSFExPSF

---
 photutils/psf/image_models.py | 108 +++++++++++++++++++++++++++++++++-
 1 file changed, 107 insertions(+), 1 deletion(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index d37616c2b..e1f1983d3 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -331,7 +331,7 @@ def evaluate(self, x, y, flux, x_0, y_0):
 
         return evaluated_model
 
-class PSFExPSF(ImagePSF):
+class PSFExPSF(Fittable2DModel):
     """
     A fittable model for PSF output from PSFEx.
 
@@ -408,8 +408,114 @@ class PSFExPSF(ImagePSF):
                     description=('Position of a feature in the image along '
                                  'the y axis'))
 
+    def __init__(self, nddata, *, flux=flux.default, x_0=x_0.default,
+                 y_0=y_0.default, origin=None, oversampling=1,
+                fill_value=0.0, **kwargs):
+
+        # self.data, self.grid_xypos = self._define_grid(nddata)
+        # self._meta = nddata.meta.copy()  # _meta to avoid the meta descriptor
+        self.oversampling = as_pair('oversampling', oversampling,
+                                    lower_bound=(0, 1))
+        self.data = nddata
+        self.fill_value = fill_value
+        # self._interpolator = {}
+
+        super().__init__(flux, x_0, y_0)
+
+    def _validate_data(data):
+        """
+        Validate the PSFEx model.
+        """
+        pass
+
+    def _cls_info(self):
+        cls_info = [('Oversampling', tuple(self.oversampling))]
+        return cls_info
+
+    def __str__(self):
+        return self._format_str(keywords=self._cls_info())
+
+    def copy(self):
+        """
+        Return a copy of this model where only the model parameters are
+        copied.
+
+        All other copied model attributes are references to the
+        original model. This prevents copying the ePSF grid data, which
+        may contain a large array.
+
+        This method is useful if one is interested in only changing
+        the model parameters in a model copy. It is used in the PSF
+        photometry classes during model fitting.
+
+        Use the `deepcopy` method if you want to copy all of the model
+        attributes, including the ePSF grid data.
+
+        Returns
+        -------
+        result : `PSFExPSF`
+            A copy of this model with only the model parameters copied.
+        """
+        newcls = object.__new__(self.__class__)
+
+        for key, val in self.__dict__.items():
+            if key in self.param_names:  # copy only the parameter values
+                newcls.__dict__[key] = copy.copy(val)
+            else:
+                newcls.__dict__[key] = val
+
+        return newcls 
+
+    def deepcopy(self):
+        """
+        Return a deep copy of this model.
+
+        Returns
+        -------
+        result : `PSFExPSF`
+            A deep copy of this model.
+        """
+        return copy.deepcopy(self)
+    
+    @property
+    def origin(self):
+        """
+        A 1D `~numpy.ndarray` (x, y) pixel coordinates within the
+        model's 2D image of the origin of the coordinate system.
+
+        The reference ``origin`` pixel will be placed at the model
+        ``x_0`` and ``y_0`` coordinates in the output coordinate system
+        on which the model is evaluated.
+
+        Most typically, the input PSF should be centered in the input
+        image, and thus the origin should be set to the central pixel of
+        the ``data`` array.
+
+        If the origin is set to `None`, then the origin will be set to
+        the center of the ``data`` array (``(npix - 1) / 2.0``).
+        """
+        return self._origin
+
+    @origin.setter
+    def origin(self, origin):
+        if origin is None:
+            origin = (np.array(self.data.shape[1:]) - 1.0) / 2.0
+        else:
+            origin = np.asarray(origin)
+            if origin.ndim != 1 or len(origin) != 2:
+                raise ValueError('origin must be 1D and have 2-elements')
+            if not np.all(np.isfinite(origin)):
+                raise ValueError('All elements of origin must be finite')
+        self._origin = origin
+
     
 
+
+
+
+
+
+
 @deprecated('2.0.0', alternative='`ImagePSF`')
 class FittableImageModel(Fittable2DModel):
     r"""

From cfad670ddbe890ca240b0e3911a67af990850c6b Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sat, 5 Oct 2024 12:08:17 +0530
Subject: [PATCH 03/14] PSFExPSF update 2

---
 photutils/psf/image_models.py | 181 +++++++++++++++++++++++++++++++++-
 1 file changed, 180 insertions(+), 1 deletion(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index b17281306..f777cea30 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -516,11 +516,190 @@ def origin(self, origin):
                 raise ValueError('All elements of origin must be finite')
         self._origin = origin
 
-    
+    @property
+    def oversampling(self):
+        """
+        The integer oversampling factor(s) of the input ePSF images.
+
+        If ``oversampling`` is a scalar then it will be used for both
+        axes. If ``oversampling`` has two elements, they must be in
+        ``(y, x)`` order.
+        """
+        return self._oversampling
+
+    @oversampling.setter
+    def oversampling(self, value):
+        """
+        Set the oversampling factor(s) of the input ePSF images.
 
+        Parameters
+        ----------
+        value : int or tuple of int
+            The integer oversampling factor(s) of the input ePSF images.
+            If ``oversampling`` is a scalar then it will be used for both
+            axes. If ``oversampling`` has two elements, they must be in
+            ``(y, x)`` order.
+        """
+        self._oversampling = as_pair('oversampling', value, lower_bound=(0, 1))
+
+    def _calc_bounding_box(self):
+        """
+        Set a bounding box defining the limits of the model.
+
+        Returns
+        -------
+        bbox : tuple
+            A bounding box defining the ((y_min, y_max), (x_min, x_max))
+            limits of the model.
+        """
+        dy, dx = np.array(self.data.shape[1:]) / 2 / self.oversampling
+        return ((self.y_0 - dy, self.y_0 + dy), (self.x_0 - dx, self.x_0 + dx))
 
+    @property
+    def bounding_box(self):
+        """
+        The bounding box of the model.
+
+        Examples
+        --------
+        >>> from itertools import product
+        >>> import numpy as np
+        >>> from astropy.nddata import NDData
+        >>> from photutils.psf import GaussianPSF, GriddedPSFModel
+        >>> psfs = []
+        >>> yy, xx = np.mgrid[0:101, 0:101]
+        >>> for i in range(16):
+        ...     theta = np.deg2rad(i * 10.0)
+        ...     gmodel = GaussianPSF(flux=1, x_0=50, y_0=50, x_fwhm=10,
+        ...                          y_fwhm=5, theta=theta)
+        ...     psfs.append(gmodel(xx, yy))
+        >>> xgrid = [0, 40, 160, 200]
+        >>> ygrid = [0, 60, 140, 200]
+        >>> meta = {}
+        >>> meta['grid_xypos'] = list(product(xgrid, ygrid))
+        >>> meta['oversampling'] = 4
+        >>> nddata = NDData(psfs, meta=meta)
+        >>> model = GriddedPSFModel(nddata, flux=1, x_0=0, y_0=0)
+        >>> model.bounding_box  # doctest: +FLOAT_CMP
+        ModelBoundingBox(
+            intervals={
+                x: Interval(lower=-12.625, upper=12.625)
+                y: Interval(lower=-12.625, upper=12.625)
+            }
+            model=GriddedPSFModel(inputs=('x', 'y'))
+            order='C'
+        )
+        """
+        return self._calc_bounding_box()
+
+
+    def _calc_interpolator(self, grid_idx):
+        """
+        Calculate the `~scipy.interpolate.RectBivariateSpline`
+        interpolator for an input ePSF image at the given reference (x,
+        y) position.
+
+        The resulting interpolator is cached in the `_interpolator`
+        dictionary for reuse.
+        """
+        xypos = tuple(self.grid_xypos[grid_idx])
+        if xypos in self._interpolator:
+            return self._interpolator[xypos]
+
+        # RectBivariateSpline expects the data to be in (x, y) axis order
+        data = self.data[grid_idx]
+        interp = RectBivariateSpline(*self._interp_xyidx, data.T, kx=3, ky=3,
+                                     s=0)
+        self._interpolator[xypos] = interp
+        return interp
 
 
+    def _calc_model_values(self, x_0, y_0, xi, yi):
+        """
+        Calculate the ePSF model at a given (x_0, y_0) model coordinate
+        and the input (xi, yi) coordinate.
+
+        Parameters
+        ----------
+        x_0, y_0 : float
+            The (x, y) position of the model.
+
+        xi, yi : float
+            The input (x, y) coordinates at which the model is
+            evaluated.
+
+        Returns
+        -------
+        result : float or `~numpy.ndarray`
+            The interpolated ePSF model at the input (x_0, y_0)
+            coordinate.
+        """
+        grid_idx, grid_xy = self._find_bounding_points(x_0, y_0)
+        interpolators = np.array([self._calc_interpolator(gidx)
+                                  for gidx in grid_idx])
+        weights = self._calc_bilinear_weights(x_0, y_0, grid_xy)
+
+        idx = np.where(weights != 0)
+        interpolators = interpolators[idx]
+        weights = weights[idx]
+
+        result = 0
+        for interp, weight in zip(interpolators, weights, strict=True):
+            result += interp(xi, yi, grid=False) * weight
+
+        return result
+
+    def evaluate(self, x, y, flux, x_0, y_0):
+        """
+        Calculate the ePSF model at the input coordinates for the given
+        model parameters.
+
+        Parameters
+        ----------
+        x, y : float or `~numpy.ndarray`
+            The x and y positions at which to evaluate the model.
+
+        flux : float
+            The flux scaling factor for the model.
+
+        x_0, y_0 : float
+            The (x, y) position of the model.
+
+        Returns
+        -------
+        evaluated_model : `~numpy.ndarray`
+            The evaluated model.
+        """
+        if x.ndim > 2:
+            raise ValueError('x and y must be 1D or 2D.')
+
+        # the base Model.__call__() method converts scalar inputs to
+        # size-1 arrays before calling evaluate(), but we need scalar
+        # values for the interpolator
+        if not np.isscalar(x_0):
+            x_0 = x_0[0]
+        if not np.isscalar(y_0):
+            y_0 = y_0[0]
+
+        # now evaluate the ePSF at the (x_0, y_0) subpixel position on
+        # the input (x, y) values
+        xi = self.oversampling[1] * (np.asarray(x, dtype=float) - x_0)
+        yi = self.oversampling[0] * (np.asarray(y, dtype=float) - y_0)
+        xi += self.origin[0]
+        yi += self.origin[1]
+
+        evaluated_model = flux * self._calc_model_values(x_0, y_0, xi, yi)
+
+        if self.fill_value is not None:
+            # set pixels that are outside the input pixel grid to the
+            # fill_value to avoid extrapolation; these bounds match the
+            # RegularGridInterpolator bounds
+            ny, nx = self.data.shape[1:]
+            invalid = (xi < 0) | (xi > nx - 1) | (yi < 0) | (yi > ny - 1)
+            evaluated_model[invalid] = self.fill_value
+
+        return evaluated_model
+
 
 
 

From ab8aa2457e871eb2d713f72d30efd84abb857c40 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sat, 5 Oct 2024 18:01:44 +0530
Subject: [PATCH 04/14] PSFExPSF update 3

---
 photutils/psf/image_models.py | 85 ++++++++++++++++++-----------------
 1 file changed, 43 insertions(+), 42 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index f777cea30..bb422df3f 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -15,6 +15,12 @@
 from photutils.aperture import CircularAperture
 from photutils.utils._parameters import as_pair
 
+try:
+    from numba import njit
+except ImportError as e:
+    def njit(func):
+        return func
+
 __all__ = ['ImagePSF', 'FittableImageModel', 'EPSFModel']
 
 
@@ -339,7 +345,7 @@ def evaluate(self, x, y, flux, x_0, y_0):
 
         return evaluated_model
 
-class PSFExPSF(Fittable2DModel):
+class PSFExVariablePSF(Fittable2DModel):
     """
     A fittable model for PSF output from PSFEx.
 
@@ -355,7 +361,7 @@ class PSFExPSF(Fittable2DModel):
 
     Parameters
     ----------
-    data : 3D `~numpy.ndarray`
+    nddata : 3D `~numpy.ndarray`
         Array containing the grid of reference PSF arrays. The length of 
         the x and y axes must both be at least 4. All elements of the 
         input image data must be finite. By default, the PSF peak is 
@@ -416,15 +422,23 @@ class PSFExPSF(Fittable2DModel):
                     description=('Position of a feature in the image along '
                                  'the y axis'))
 
-    def __init__(self, nddata, *, flux=flux.default, x_0=x_0.default,
+    def __init__(self, nddata, meta, *, flux=flux.default, x_0=x_0.default,
                  y_0=y_0.default, origin=None, oversampling=1,
                 fill_value=0.0, **kwargs):
 
         # self.data, self.grid_xypos = self._define_grid(nddata)
         # self._meta = nddata.meta.copy()  # _meta to avoid the meta descriptor
-        self.oversampling = as_pair('oversampling', oversampling,
-                                    lower_bound=(0, 1))
+
         self.data = nddata
+        self.psf_shape = self.data[0].shape
+        self._meta = meta
+        self.vardeg = meta['POLDEG1']
+        self.xzero = meta['POLZERO1']
+        self.yzero = meta['POLZERO2']
+        self.xscale = meta['POLSCAL1']
+        self.yscale = meta['POLSCAL2']
+        self.oversampling = as_pair('oversampling', meta['PSF_SAMP'],
+                                    lower_bound=(0, 1))
         self.fill_value = fill_value
         # self._interpolator = {}
 
@@ -433,6 +447,11 @@ def __init__(self, nddata, *, flux=flux.default, x_0=x_0.default,
     def _validate_data(data):
         """
         Validate the PSFEx model.
+
+        1. Oversampling should be interger
+        2. first independent variable should be X_IMAGE
+        3. Second independent variable should be Y_IMAGE
+        4. Both variables should be in Group 1.
         """
         pass
 
@@ -552,7 +571,7 @@ def _calc_bounding_box(self):
             A bounding box defining the ((y_min, y_max), (x_min, x_max))
             limits of the model.
         """
-        dy, dx = np.array(self.data.shape[1:]) / 2 / self.oversampling
+        dy, dx = np.array(self.psf_shape) / 2 / self.oversampling
         return ((self.y_0 - dy, self.y_0 + dy), (self.x_0 - dx, self.x_0 + dx))
 
     @property
@@ -592,49 +611,31 @@ def bounding_box(self):
         """
         return self._calc_bounding_box()
 
+    @njit
+    def _calc_poly_coeffs(x, y, vardeg):
+        coeffs = []
+        for i in range(vardeg+1):
+            for j in range(vardeg-i+1):
+                print(f'x^{j}y^{i}')
+                coeffs.append(x**j * y**i)
+        return coeffs
 
-    def _calc_interpolator(self, grid_idx):
-        """
-        Calculate the `~scipy.interpolate.RectBivariateSpline`
-        interpolator for an input ePSF image at the given reference (x,
-        y) position.
-
-        The resulting interpolator is cached in the `_interpolator`
-        dictionary for reuse.
-        """
-        xypos = tuple(self.grid_xypos[grid_idx])
-        if xypos in self._interpolator:
-            return self._interpolator[xypos]
-
-        # RectBivariateSpline expects the data to be in (x, y) axis order
-        data = self.data[grid_idx]
-        interp = RectBivariateSpline(*self._interp_xyidx, data.T, kx=3, ky=3,
-                                     s=0)
-        self._interpolator[xypos] = interp
-        return interp
 
+    def _calc_image_weights(self, x, y):
+        xi = (x - self.xzero) / self.xscale
+        yi = (y - self.yzero) / self.yscale
+        return self._calc_poly_coeffs(xi, yi, self.vardeg)
 
-    def _calc_model_values(self, x_0, y_0, xi, yi):
+    def _calc_model_values(self):
+        """
+        Docs to be updated
         """
-        Calculate the ePSF model at a given (x_0, y_0) model coordinate
-        and the input (xi, yi) coordinate.
 
-        Parameters
-        ----------
-        x_0, y_0 : float
-            The (x, y) position of the model.
 
-        xi, yi : float
-            The input (x, y) coordinates at which the model is
-            evaluated.
+        psf = sum(a*b for a, b in zip(coeffs, psf_model))
+        psf_models.append(psf)
+
 
-        Returns
-        -------
-        result : float or `~numpy.ndarray`
-            The interpolated ePSF model at the input (x_0, y_0)
-            coordinate.
-        """
-        grid_idx, grid_xy = self._find_bounding_points(x_0, y_0)
         interpolators = np.array([self._calc_interpolator(gidx)
                                   for gidx in grid_idx])
         weights = self._calc_bilinear_weights(x_0, y_0, grid_xy)

From 9f80b7199e77607036738c94bf9f77a1e350cbf3 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sat, 5 Oct 2024 18:59:30 +0530
Subject: [PATCH 05/14] Complete rough draft if PSFExPSF

---
 photutils/psf/image_models.py | 54 ++++++++++++++++-------------------
 1 file changed, 24 insertions(+), 30 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index bb422df3f..66c8ad0c3 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -620,35 +620,41 @@ def _calc_poly_coeffs(x, y, vardeg):
                 coeffs.append(x**j * y**i)
         return coeffs
 
-
     def _calc_image_weights(self, x, y):
         xi = (x - self.xzero) / self.xscale
         yi = (y - self.yzero) / self.yscale
         return self._calc_poly_coeffs(xi, yi, self.vardeg)
 
-    def _calc_model_values(self):
-        """
-        Docs to be updated
+    def _calc_model_values(self, x_0, y_0, xi, yi):
         """
+        Calculate the ePSF model at a given (x_0, y_0) model coordinate
+        and the input (xi, yi) coordinate.
 
+        Parameters
+        ----------
+        x_0, y_0 : float
+            The (x, y) position of the model.
 
-        psf = sum(a*b for a, b in zip(coeffs, psf_model))
-        psf_models.append(psf)
-
+        xi, yi : float
+            The input (x, y) coordinates at which the model is
+            evaluated.
 
-        interpolators = np.array([self._calc_interpolator(gidx)
-                                  for gidx in grid_idx])
-        weights = self._calc_bilinear_weights(x_0, y_0, grid_xy)
+        Returns
+        -------
+        result : float or `~numpy.ndarray`
+            The interpolated ePSF model at the input (x_0, y_0)
+            coordinate.
+        """
 
-        idx = np.where(weights != 0)
-        interpolators = interpolators[idx]
-        weights = weights[idx]
+        weights = self._calc_image_weights(x_0, y_0)
+        psf = sum(a*b for a, b in zip(weights, self.data))
 
-        result = 0
-        for interp, weight in zip(interpolators, weights, strict=True):
-            result += interp(xi, yi, grid=False) * weight
+        x = np.arange(self.psf_shape[1])
+        y = np.arange(self.psf_shape[0])
+        # RectBivariateSpline expects the data to be in (x, y) axis order
+        interpolator = RectBivariateSpline(x, y, psf.T, kx=3, ky=3, s=0)
 
-        return result
+        return interpolator(xi, yi, grid=False)
 
     def evaluate(self, x, y, flux, x_0, y_0):
         """
@@ -671,19 +677,7 @@ def evaluate(self, x, y, flux, x_0, y_0):
         evaluated_model : `~numpy.ndarray`
             The evaluated model.
         """
-        if x.ndim > 2:
-            raise ValueError('x and y must be 1D or 2D.')
-
-        # the base Model.__call__() method converts scalar inputs to
-        # size-1 arrays before calling evaluate(), but we need scalar
-        # values for the interpolator
-        if not np.isscalar(x_0):
-            x_0 = x_0[0]
-        if not np.isscalar(y_0):
-            y_0 = y_0[0]
 
-        # now evaluate the ePSF at the (x_0, y_0) subpixel position on
-        # the input (x, y) values
         xi = self.oversampling[1] * (np.asarray(x, dtype=float) - x_0)
         yi = self.oversampling[0] * (np.asarray(y, dtype=float) - y_0)
         xi += self.origin[0]
@@ -695,7 +689,7 @@ def evaluate(self, x, y, flux, x_0, y_0):
             # set pixels that are outside the input pixel grid to the
             # fill_value to avoid extrapolation; these bounds match the
             # RegularGridInterpolator bounds
-            ny, nx = self.data.shape[1:]
+            ny, nx = self.psf_shape
             invalid = (xi < 0) | (xi > nx - 1) | (yi < 0) | (yi > ny - 1)
             evaluated_model[invalid] = self.fill_value
 

From d95791d9888f1225f7742ff416e8887dd25c2938 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 01:33:39 +0530
Subject: [PATCH 06/14] Clean up the class for initial review

---
 photutils/psf/image_models.py | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 66c8ad0c3..7eed7ab05 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -352,7 +352,7 @@ class PSFExVariablePSF(Fittable2DModel):
     The PSFEx PSF is a grid of 2D models where each 2D element represents
     the contribution of a polynomial. The polymial can be constructed from
     any FITS header parameter, but the most straightforward example would
-    be the x and y positions of the CCD image.
+    be the x and y positions of the CCD image. Currently, only variation in x and y is supported.
 
     The model has three model parameters: an image intensity scaling
     factor (``flux``) which is applied to the input image, and two
@@ -368,6 +368,9 @@ class PSFExVariablePSF(Fittable2DModel):
         assumed to be located at the center of the input image (see the ``origin``
         keyword).
 
+    meta : `~astropy.io.fits.header.Header`
+        Header associated with extension PSFEx PSF.
+
     flux : float, optional
         The total flux of the source, assuming the input image
         was properly normalized.
@@ -399,6 +402,8 @@ class PSFExVariablePSF(Fittable2DModel):
         ``oversampling`` is a scalar then it will be used for both axes.
         If ``oversampling`` has two elements, they must be in ``(y, x)``
         order.
+        PSFEx supports oversampling as a float, but photutils supports 
+        it as an interger as of now.
 
     fill_value : float, optional
         The value to use for points outside of the input pixel grid.
@@ -410,6 +415,7 @@ class PSFExVariablePSF(Fittable2DModel):
 
     See Also
     --------
+    ImagePSF : A model for a 2D image PSF.
     GriddedPSFModel : A model for a grid of ePSF models.    
 
     """
@@ -426,11 +432,8 @@ def __init__(self, nddata, meta, *, flux=flux.default, x_0=x_0.default,
                  y_0=y_0.default, origin=None, oversampling=1,
                 fill_value=0.0, **kwargs):
 
-        # self.data, self.grid_xypos = self._define_grid(nddata)
-        # self._meta = nddata.meta.copy()  # _meta to avoid the meta descriptor
-
         self.data = nddata
-        self.psf_shape = self.data[0].shape
+        self.psf_shape = self.data.shape[1:]
         self._meta = meta
         self.vardeg = meta['POLDEG1']
         self.xzero = meta['POLZERO1']
@@ -439,6 +442,7 @@ def __init__(self, nddata, meta, *, flux=flux.default, x_0=x_0.default,
         self.yscale = meta['POLSCAL2']
         self.oversampling = as_pair('oversampling', meta['PSF_SAMP'],
                                     lower_bound=(0, 1))
+        self.origin = origin
         self.fill_value = fill_value
         # self._interpolator = {}
 
@@ -480,7 +484,7 @@ def copy(self):
 
         Returns
         -------
-        result : `PSFExPSF`
+        result : `PSFExVariablePSF`
             A copy of this model with only the model parameters copied.
         """
         newcls = object.__new__(self.__class__)
@@ -499,7 +503,7 @@ def deepcopy(self):
 
         Returns
         -------
-        result : `PSFExPSF`
+        result : `PSFExVariablePSF`
             A deep copy of this model.
         """
         return copy.deepcopy(self)
@@ -526,7 +530,7 @@ def origin(self):
     @origin.setter
     def origin(self, origin):
         if origin is None:
-            origin = (np.array(self.data.shape[1:]) - 1.0) / 2.0
+            origin = (np.array(self.psf_shape) - 1.0) / 2.0
         else:
             origin = np.asarray(origin)
             if origin.ndim != 1 or len(origin) != 2:
@@ -616,7 +620,6 @@ def _calc_poly_coeffs(x, y, vardeg):
         coeffs = []
         for i in range(vardeg+1):
             for j in range(vardeg-i+1):
-                print(f'x^{j}y^{i}')
                 coeffs.append(x**j * y**i)
         return coeffs
 
@@ -696,8 +699,6 @@ def evaluate(self, x, y, flux, x_0, y_0):
         return evaluated_model
 
 
-
-
 @deprecated('2.0.0', alternative='`ImagePSF`')
 class FittableImageModel(Fittable2DModel):
     r"""

From b6c2e23710df8ea8cf53feafbbd55a3c52ef8567 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 02:26:37 +0530
Subject: [PATCH 07/14] Fix formatting

---
 photutils/psf/image_models.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 7eed7ab05..3a6c281b2 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -345,6 +345,7 @@ def evaluate(self, x, y, flux, x_0, y_0):
 
         return evaluated_model
 
+
 class PSFExVariablePSF(Fittable2DModel):
     """
     A fittable model for PSF output from PSFEx.
@@ -352,7 +353,8 @@ class PSFExVariablePSF(Fittable2DModel):
     The PSFEx PSF is a grid of 2D models where each 2D element represents
     the contribution of a polynomial. The polymial can be constructed from
     any FITS header parameter, but the most straightforward example would
-    be the x and y positions of the CCD image. Currently, only variation in x and y is supported.
+    be the x and y positions of the CCD image. Currently, only variation in x 
+    and y is supported.
 
     The model has three model parameters: an image intensity scaling
     factor (``flux``) which is applied to the input image, and two
@@ -365,8 +367,8 @@ class PSFExVariablePSF(Fittable2DModel):
         Array containing the grid of reference PSF arrays. The length of 
         the x and y axes must both be at least 4. All elements of the 
         input image data must be finite. By default, the PSF peak is 
-        assumed to be located at the center of the input image (see the ``origin``
-        keyword).
+        assumed to be located at the center of the input image (see the 
+        ``origin`` keyword).
 
     meta : `~astropy.io.fits.header.Header`
         Header associated with extension PSFEx PSF.
@@ -444,10 +446,10 @@ def __init__(self, nddata, meta, *, flux=flux.default, x_0=x_0.default,
                                     lower_bound=(0, 1))
         self.origin = origin
         self.fill_value = fill_value
-        # self._interpolator = {}
 
         super().__init__(flux, x_0, y_0)
 
+    @staticmethod
     def _validate_data(data):
         """
         Validate the PSFEx model.
@@ -460,8 +462,7 @@ def _validate_data(data):
         pass
 
     def _cls_info(self):
-        cls_info = [('Oversampling', tuple(self.oversampling))]
-        return cls_info
+        return [('Oversampling', tuple(self.oversampling))]
 
     def __str__(self):
         return self._format_str(keywords=self._cls_info())

From 71da9a7927f4bd594c9c9bea3dc9f4c940d53369 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 02:31:38 +0530
Subject: [PATCH 08/14] Fix PEP8 warnings

---
 photutils/psf/image_models.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 3a6c281b2..952f21611 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -353,7 +353,7 @@ class PSFExVariablePSF(Fittable2DModel):
     The PSFEx PSF is a grid of 2D models where each 2D element represents
     the contribution of a polynomial. The polymial can be constructed from
     any FITS header parameter, but the most straightforward example would
-    be the x and y positions of the CCD image. Currently, only variation in x 
+    be the x and y positions of the CCD image. Currently, only variation in x
     and y is supported.
 
     The model has three model parameters: an image intensity scaling
@@ -364,10 +364,10 @@ class PSFExVariablePSF(Fittable2DModel):
     Parameters
     ----------
     nddata : 3D `~numpy.ndarray`
-        Array containing the grid of reference PSF arrays. The length of 
-        the x and y axes must both be at least 4. All elements of the 
-        input image data must be finite. By default, the PSF peak is 
-        assumed to be located at the center of the input image (see the 
+        Array containing the grid of reference PSF arrays. The length of
+        the x and y axes must both be at least 4. All elements of the
+        input image data must be finite. By default, the PSF peak is
+        assumed to be located at the center of the input image (see the
         ``origin`` keyword).
 
     meta : `~astropy.io.fits.header.Header`
@@ -404,7 +404,7 @@ class PSFExVariablePSF(Fittable2DModel):
         ``oversampling`` is a scalar then it will be used for both axes.
         If ``oversampling`` has two elements, they must be in ``(y, x)``
         order.
-        PSFEx supports oversampling as a float, but photutils supports 
+        PSFEx supports oversampling as a float, but photutils supports
         it as an interger as of now.
 
     fill_value : float, optional
@@ -418,7 +418,7 @@ class PSFExVariablePSF(Fittable2DModel):
     See Also
     --------
     ImagePSF : A model for a 2D image PSF.
-    GriddedPSFModel : A model for a grid of ePSF models.    
+    GriddedPSFModel : A model for a grid of ePSF models.
 
     """
     flux = Parameter(default=1,
@@ -432,7 +432,7 @@ class PSFExVariablePSF(Fittable2DModel):
 
     def __init__(self, nddata, meta, *, flux=flux.default, x_0=x_0.default,
                  y_0=y_0.default, origin=None, oversampling=1,
-                fill_value=0.0, **kwargs):
+                 fill_value=0.0, **kwargs):
 
         self.data = nddata
         self.psf_shape = self.data.shape[1:]
@@ -496,7 +496,7 @@ def copy(self):
             else:
                 newcls.__dict__[key] = val
 
-        return newcls 
+        return newcls
 
     def deepcopy(self):
         """
@@ -508,7 +508,7 @@ def deepcopy(self):
             A deep copy of this model.
         """
         return copy.deepcopy(self)
-    
+
     @property
     def origin(self):
         """

From 17514f51c5c42ce423bb3b09875992ef7cf8edeb Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 13:01:44 +0530
Subject: [PATCH 09/14] Attempt 1 at fixing pre-commit warnings

---
 photutils/psf/image_models.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 952f21611..071b08b42 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -616,11 +616,12 @@ def bounding_box(self):
         """
         return self._calc_bounding_box()
 
+    @staticmethod
     @njit
     def _calc_poly_coeffs(x, y, vardeg):
         coeffs = []
-        for i in range(vardeg+1):
-            for j in range(vardeg-i+1):
+        for i in range(vardeg + 1):
+            for j in range(vardeg - i + 1):
                 coeffs.append(x**j * y**i)
         return coeffs
 
@@ -651,7 +652,7 @@ def _calc_model_values(self, x_0, y_0, xi, yi):
         """
 
         weights = self._calc_image_weights(x_0, y_0)
-        psf = sum(a*b for a, b in zip(weights, self.data))
+        psf = sum(a * b for a, b in zip(weights, self.data))
 
         x = np.arange(self.psf_shape[1])
         y = np.arange(self.psf_shape[0])

From 4348307155274f0a8f000b89976862efcddb8967 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 13:08:28 +0530
Subject: [PATCH 10/14] Attempt 2 at fixing pre-commit warnings

---
 photutils/psf/image_models.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 071b08b42..2324dcb32 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -419,7 +419,6 @@ class PSFExVariablePSF(Fittable2DModel):
     --------
     ImagePSF : A model for a 2D image PSF.
     GriddedPSFModel : A model for a grid of ePSF models.
-
     """
     flux = Parameter(default=1,
                      description='Intensity scaling factor of the image.')

From c4d8259cf84532187449370952754346de6f4eff Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 16:47:46 +0530
Subject: [PATCH 11/14] Attempt 3 at fixing pre-commit warnings

---
 photutils/psf/image_models.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 2324dcb32..8e59ab45c 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -618,11 +618,9 @@ def bounding_box(self):
     @staticmethod
     @njit
     def _calc_poly_coeffs(x, y, vardeg):
-        coeffs = []
-        for i in range(vardeg + 1):
-            for j in range(vardeg - i + 1):
-                coeffs.append(x**j * y**i)
-        return coeffs
+        return [x**j * y**i 
+                for i in range(vardeg+1) 
+                for j in range(vardeg-i+1)]
 
     def _calc_image_weights(self, x, y):
         xi = (x - self.xzero) / self.xscale

From 1e3d1af5b7e5a0eaf83c54f85445ba9b129effa8 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 16:53:28 +0530
Subject: [PATCH 12/14] Attempt 4 at fixing pre-commit warnings

---
 photutils/psf/image_models.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 8e59ab45c..7b5691040 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -618,8 +618,8 @@ def bounding_box(self):
     @staticmethod
     @njit
     def _calc_poly_coeffs(x, y, vardeg):
-        return [x**j * y**i 
-                for i in range(vardeg+1) 
+        return [x**j * y**i
+                for i in range(vardeg+1)
                 for j in range(vardeg-i+1)]
 
     def _calc_image_weights(self, x, y):

From 1d968769712cd9bbbcbb0de4bb675c045ee3b653 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 16:54:48 +0530
Subject: [PATCH 13/14] Attempt 5 at fixing pre-commit warnings

---
 photutils/psf/image_models.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index 7b5691040..d7216afdd 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -619,8 +619,8 @@ def bounding_box(self):
     @njit
     def _calc_poly_coeffs(x, y, vardeg):
         return [x**j * y**i
-                for i in range(vardeg+1)
-                for j in range(vardeg-i+1)]
+                for i in range(vardeg + 1)
+                for j in range(vardeg - i + 1)]
 
     def _calc_image_weights(self, x, y):
         xi = (x - self.xzero) / self.xscale

From 53d675fc9c8e1d71c5ab2c96db016a140366fb78 Mon Sep 17 00:00:00 2001
From: Naveen Dukiya <97835976+navii98@users.noreply.github.com>
Date: Sun, 6 Oct 2024 20:39:21 +0530
Subject: [PATCH 14/14] Implement validate_data method

---
 photutils/psf/image_models.py | 51 ++++++++++++++++++++++++++---------
 1 file changed, 38 insertions(+), 13 deletions(-)

diff --git a/photutils/psf/image_models.py b/photutils/psf/image_models.py
index d7216afdd..3a8a3c8e0 100644
--- a/photutils/psf/image_models.py
+++ b/photutils/psf/image_models.py
@@ -354,7 +354,8 @@ class PSFExVariablePSF(Fittable2DModel):
     the contribution of a polynomial. The polymial can be constructed from
     any FITS header parameter, but the most straightforward example would
     be the x and y positions of the CCD image. Currently, only variation in x
-    and y is supported.
+    and y is supported. e.g. for a variation degree of 2, the output PSF will
+    be a grid of 6 images corresponding to 1, x, x^2, y, xy, and y^2.
 
     The model has three model parameters: an image intensity scaling
     factor (``flux``) which is applied to the input image, and two
@@ -433,6 +434,7 @@ def __init__(self, nddata, meta, *, flux=flux.default, x_0=x_0.default,
                  y_0=y_0.default, origin=None, oversampling=1,
                  fill_value=0.0, **kwargs):
 
+        self._validate_data(data, meta)
         self.data = nddata
         self.psf_shape = self.data.shape[1:]
         self._meta = meta
@@ -441,27 +443,46 @@ def __init__(self, nddata, meta, *, flux=flux.default, x_0=x_0.default,
         self.yzero = meta['POLZERO2']
         self.xscale = meta['POLSCAL1']
         self.yscale = meta['POLSCAL2']
-        self.oversampling = as_pair('oversampling', meta['PSF_SAMP'],
+        self.oversampling = as_pair('oversampling', int(meta['PSF_SAMP']),
                                     lower_bound=(0, 1))
         self.origin = origin
         self.fill_value = fill_value
 
-        super().__init__(flux, x_0, y_0)
+        super().__init__(flux, x_0, y_0, **kwargs)
 
     @staticmethod
-    def _validate_data(data):
-        """
-        Validate the PSFEx model.
+    def _validate_data(data, meta):
+        if not isinstance(data, np.ndarray):
+            raise TypeError('Input data must be a 3D numpy array.')
 
-        1. Oversampling should be interger
-        2. first independent variable should be X_IMAGE
-        3. Second independent variable should be Y_IMAGE
-        4. Both variables should be in Group 1.
-        """
-        pass
+        if data.ndim != 3:
+            raise ValueError('Input data must be a 3D numpy array.')
+
+        if not np.all(np.isfinite(data)):
+            raise ValueError('All elements of input data must be finite.')
+
+        # this is required by RectBivariateSpline for kx=3, ky=3
+        if np.any(np.array(data.shape[1:]) < 4):
+            raise ValueError('The length of the x and y axes must both be at '
+                             'least 4.')
+        
+        if not meta['PSF_SAMP'].is_integer():
+            raise TypeError('Oversampling is not supported as a flaot.')
+
+        if (('POLNAXIS' in meta.keys() and meta['POLNAXIS'] > 2) or
+           ('POLNGRP' in meta.keys() and meta['POLNGRP'] > 1) or
+           ('POLNAME1' in meta.keys() and meta['POLNAME1'] != 'X_IMAGE') or
+           ('POLNAME2' in meta.keys() and meta['POLNAME2'] != 'Y_IMAGE')
+        ):
+            raise NotImplementedError('Only one variability group is'
+                                      'supported with upto two variables'
+                                      'that should be "X_IMAGE" and "Y_IMAGE",'
+                                      'in that order.')
 
     def _cls_info(self):
-        return [('Oversampling', tuple(self.oversampling))]
+        return [('Oversampling', tuple(self.oversampling)),
+                ('PSF Shape (oversampled pixels)', self.psf_shape),
+                ('Variation degree', self.vardeg)]
 
     def __str__(self):
         return self._format_str(keywords=self._cls_info())
@@ -623,6 +644,10 @@ def _calc_poly_coeffs(x, y, vardeg):
                 for j in range(vardeg - i + 1)]
 
     def _calc_image_weights(self, x, y):
+        """
+        Calculates the weights for each 2D array of PSF grid based on the
+        location of the source in output grid (x, y).
+        """
         xi = (x - self.xzero) / self.xscale
         yi = (y - self.yzero) / self.yscale
         return self._calc_poly_coeffs(xi, yi, self.vardeg)