From f448f03f249deec722f7f1e075e32992a43ea18d Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Wed, 3 May 2023 13:55:49 -0400
Subject: [PATCH 01/12] Initial implementation of SumPSF

---
 piff/__init__.py     |   3 +
 piff/psf.py          |  42 ++++--
 piff/simplepsf.py    |   5 -
 piff/star.py         |   2 +-
 piff/sumpsf.py       | 312 +++++++++++++++++++++++++++++++++++++++++++
 tests/test_sumpsf.py | 147 ++++++++++++++++++++
 6 files changed, 491 insertions(+), 20 deletions(-)
 create mode 100644 piff/sumpsf.py
 create mode 100644 tests/test_sumpsf.py

diff --git a/piff/__init__.py b/piff/__init__.py
index 16f3fb90..f8c49de1 100644
--- a/piff/__init__.py
+++ b/piff/__init__.py
@@ -111,6 +111,9 @@
 # Optics
 from .optical_model import Optical, optical_templates
 
+# Composite PSFs
+from .sumpsf import SumPSF
+
 # Leave these in their own namespaces
 from . import util
 from . import des
diff --git a/piff/psf.py b/piff/psf.py
index 8a6d899e..6c81b149 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -101,6 +101,11 @@ def set_num(self, num):
         # But this is the minimum action that all subclasses need to do.
         self._num = num
 
+    @property
+    def num_components(self):
+        # Subclasses for which this is not true can overwrite this
+        return 1
+
     @classmethod
     def __init_subclass__(cls):
         # Classes that don't want to register a type name can either not define _type_name
@@ -694,7 +699,12 @@ def _drawStar(self, star, center=None):
         raise NotImplementedError("Derived classes must define the _drawStar function")
 
     def _getProfile(self, star):
-        raise NotImplementedError("Derived classes must define the _getProfile function")
+        prof, method = self._getRawProfile(star)
+        prof = prof.shift(star.fit.center) * star.fit.flux
+        return prof, method
+
+    def _getRawProfile(self, star):
+        raise NotImplementedError("Derived classes must define the _getRawProfile function")
 
     def write(self, file_name, logger=None):
         """Write a PSF object to a file.
@@ -723,10 +733,12 @@ def _write(self, fits, extname, logger=None):
         psf_type = self._type_name
         write_kwargs(fits, extname, dict(self.kwargs, type=psf_type, piff_version=piff_version))
         logger.info("Wrote the basic PSF information to extname %s", extname)
-        Star.write(self.stars, fits, extname=extname + '_stars')
-        logger.info("Wrote the PSF stars to extname %s", extname + '_stars')
-        self.writeWCS(fits, extname=extname + '_wcs', logger=logger)
-        logger.info("Wrote the PSF WCS to extname %s", extname + '_wcs')
+        if hasattr(self, 'stars'):
+            Star.write(self.stars, fits, extname=extname + '_stars')
+            logger.info("Wrote the PSF stars to extname %s", extname + '_stars')
+        if hasattr(self, 'wcs'):
+            self.writeWCS(fits, extname=extname + '_wcs', logger=logger)
+            logger.info("Wrote the PSF WCS to extname %s", extname + '_wcs')
         self._finish_write(fits, extname=extname, logger=logger)
 
     @classmethod
@@ -773,17 +785,19 @@ def _read(cls, fits, extname, logger):
             raise ValueError("psf type %s is not a valid Piff PSF"%psf_type)
         psf_cls = PSF.valid_psf_types[psf_type]
 
-        # Read the stars, wcs, pointing values
-        stars = Star.read(fits, extname + '_stars')
-        logger.debug("stars = %s",stars)
-        wcs, pointing = cls.readWCS(fits, extname + '_wcs', logger=logger)
-        logger.debug("wcs = %s, pointing = %s",wcs,pointing)
-
         # Make the PSF instance
         psf = psf_cls(**kwargs)
-        psf.stars = stars
-        psf.wcs = wcs
-        psf.pointing = pointing
+
+        # Read the stars, wcs, pointing values
+        if extname + '_stars' in fits:
+            stars = Star.read(fits, extname + '_stars')
+            logger.debug("stars = %s",stars)
+            psf.stars = stars
+        if extname + '_wcs' in fits:
+            wcs, pointing = cls.readWCS(fits, extname + '_wcs', logger=logger)
+            logger.debug("wcs = %s, pointing = %s",wcs,pointing)
+            psf.wcs = wcs
+            psf.pointing = pointing
 
         # Just in case the class needs to do something else at the end.
         psf._finish_read(fits, extname, logger)
diff --git a/piff/simplepsf.py b/piff/simplepsf.py
index 507cc04e..87ced184 100644
--- a/piff/simplepsf.py
+++ b/piff/simplepsf.py
@@ -228,11 +228,6 @@ def interpolateStar(self, star):
     def _drawStar(self, star, center=None):
         return self.model.draw(star, center=center)
 
-    def _getProfile(self, star):
-        prof, method = self._getRawProfile(star)
-        prof = prof.shift(star.fit.center) * star.fit.flux
-        return prof, method
-
     def _getRawProfile(self, star):
         return self.model.getProfile(star.fit.get_params(self._num)), self.model._method
 
diff --git a/piff/star.py b/piff/star.py
index 2196c612..0d070b48 100644
--- a/piff/star.py
+++ b/piff/star.py
@@ -384,7 +384,7 @@ def write(self, stars, fits, extname):
 
             # params might need to be flattened
             if stars[0].fit.params_lens is not None:
-                header = {'PARAMS_LENS' : stars[0].fit.params_lens}
+                header = {'PARAMS_LENS' : str(stars[0].fit.params_lens)}
                 params = [ StarFit.flatten_params(p) for p in params ]
 
             dtypes.append( ('params', float, len(params[0])) )
diff --git a/piff/sumpsf.py b/piff/sumpsf.py
new file mode 100644
index 00000000..2bd9f4a5
--- /dev/null
+++ b/piff/sumpsf.py
@@ -0,0 +1,312 @@
+# Copyright (c) 2016 by Mike Jarvis and the other collaborators on GitHub at
+# https://github.com/rmjarvis/Piff  All rights reserved.
+#
+# Piff is free software: Redistribution and use in source and binary forms
+# with or without modification, are permitted provided that the following
+# conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the disclaimer given in the accompanying LICENSE
+#    file.
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the disclaimer given in the documentation
+#    and/or other materials provided with the distribution.
+
+"""
+.. module:: psf
+"""
+
+import numpy as np
+import galsim
+
+from .psf import PSF
+from .util import write_kwargs, read_kwargs
+from .star import Star, StarFit
+from .outliers import Outliers
+
+class SumPSF(PSF):
+    """A PSF class that is the Sum of two or more other PSF types.
+
+    A SumPSF is built from an ordered list of other PSFs.
+
+    When fitting the Sum, the default pattern is that all components after the first one
+    are initialized to zero, and the first component is fit as usual, but just using a
+    single iteration of the fit.  Then the residuals of this model are fit using the second
+    component, and so on.  Once all components are fit, outliers may be rejected, and then
+    the process is iterated.
+
+    This pattern can be tweaked somewhat using the initialization options available to
+    PSF models.  If a component should be initialized to something other than a zero model,
+    then one should explicitly set it.
+
+    :param components:  A list of PSF instances defining the components to be summed.
+    :param outliers:    Optionally, an Outliers instance used to remove outliers.
+                        [default: None]
+    :param chisq_thresh: Change in reduced chisq at which iteration will terminate.
+                        [default: 0.1]
+    :param max_iter:    Maximum number of iterations to try. [default: 30]
+    """
+    _type_name = 'Sum'
+
+    def __init__(self, components, outliers=None, chisq_thresh=0.1, max_iter=30):
+        self.components = components
+        self.outliers = outliers
+        self.chisq_thresh = chisq_thresh
+        self.max_iter = max_iter
+        self.kwargs = {
+            # If components is a list, mark the number of components here for I/O purposed.
+            # But if it's just a number, leave it alone.
+            'components': len(components) if isinstance(components, list) else components,
+            'outliers': 0,
+            'chisq_thresh': self.chisq_thresh,
+            'max_iter': self.max_iter,
+        }
+        self.chisq = 0.
+        self.last_delta_chisq = 0.
+        self.dof = 0
+        self.nremoved = 0
+        self.niter = 0
+        self.set_num(None)
+
+    def set_num(self, num):
+        """If there are multiple components involved in the fit, set the number to use
+        for this model.
+        """
+        if isinstance(self.components, list):
+            # Then building components has been completed.
+
+            # This array keeps track of which num to use for each component.
+            self._nums = np.empty(len(self.components), dtype=int)
+            self._num_components = 0  # It might not be len(self.components) if some of them are
+                                      # in turn composite types.  Figure it out below.
+
+            k1 = 0 if num is None else num
+            for k, comp in enumerate(self.components):
+                self._nums[k] = k1
+                comp.set_num(k1)
+                k1 += comp.num_components
+                self._num_components += comp.num_components
+            self._num = self._nums[0]
+        else:
+            self._num = None
+
+        # else components are not yet built. This will be called again when that's done.
+
+    @property
+    def num_components(self):
+        return self._num_components
+
+    @property
+    def interp_property_names(self):
+        names = set()
+        for c in self.components:
+            names.update(c.interp_property_names)
+        return names
+
+    @classmethod
+    def parseKwargs(cls, config_psf, logger):
+        """Parse the psf field of a configuration dict and return the kwargs to use for
+        initializing an instance of the class.
+
+        :param config_psf:      The psf field of the configuration dict, config['psf']
+        :param logger:          A logger object for logging debug info. [default: None]
+
+        :returns: a kwargs dict to pass to the initializer
+        """
+        from .outliers import Outliers
+
+        kwargs = {}
+        kwargs.update(config_psf)
+        kwargs.pop('type',None)
+
+        if 'components' not in kwargs:
+            raise ValueError("components field is required in psf field for type=Sum")
+
+        # make components
+        components = kwargs.pop('components')
+        kwargs['components'] = []
+        for comp in components:
+            kwargs['components'].append(PSF.process(comp, logger=logger))
+
+        if 'outliers' in kwargs:
+            outliers = Outliers.process(kwargs.pop('outliers'), logger=logger)
+            kwargs['outliers'] = outliers
+
+        return kwargs
+
+    def setup_params(self, stars):
+        """Make sure the stars have the right shape params object
+        """
+        new_stars = []
+        for star in stars:
+            if star.fit.params is None:
+                fit = star.fit.withNew(params=[None] * self.num_components,
+                                       params_var=[None] * self.num_components)
+                star = Star(star.data, fit)
+            else:
+                assert len(star.fit.params) > self._nums[-1]
+            new_stars.append(star)
+        return new_stars
+
+    def initialize_params(self, stars, logger, default_init=None):
+        nremoved = 0
+
+        logger.debug("Initializing components of SumPSF")
+
+        # First make sure the params are the right shape.
+        stars = self.setup_params(stars)
+
+        # Now initialize all the components
+        for k, comp in enumerate(self.components):
+            stars, nremoved1 = comp.initialize_params(stars, logger, default_init=default_init)
+            nremoved += nremoved1
+            # After the first one, set default_init to 'zero'
+            default_init='zero'
+
+        # If any components are degenerate, mark the sum as degenerate.
+        self.degenerate_points = any([comp.degenerate_points for comp in self.components])
+
+        return stars, nremoved
+
+    def single_iteration(self, stars, logger, convert_func):
+        nremoved = 0  # For this iteration
+
+        # Draw the current models for each component
+        current_models = []
+        for k, comp in enumerate(self.components):
+            current_models.append(comp.drawStarList(stars))
+
+        # Fit each component in order
+        for k, comp in enumerate(self.components):
+            # Update stars to subtract current fit from other components.
+            modified_stars = []
+            for i, star in enumerate(stars):
+                new_image = star.image.copy()
+                for kk in range(len(self.components)):
+                    if kk != k:
+                        new_image -= current_models[kk][i].image
+                modified_star = Star(star.data.withNew(image=new_image), star.fit)
+                modified_stars.append(modified_star)
+
+            # Fit this component
+            new_stars, nremoved1 = comp.single_iteration(modified_stars, logger, convert_func)
+
+            # The fit components in new_stars are what we want, but the data parts are
+            # corrupted by subtracting the other components, which we don't want.
+            # Just copy over the fits.
+            stars = [Star(s1.data, s2.fit) for s1,s2 in zip(stars, new_stars)]
+            nremoved += nremoved1
+
+            # Update the current models for later components
+            stars = comp.interpolateStarList(stars)
+            current_models[k] = comp.drawStarList(stars)
+
+        return stars, nremoved
+
+    @property
+    def fit_center(self):
+        """Whether to fit the center of the star in reflux.
+
+        This is generally set in the model specifications.
+        If all component models includes a shift, then this is False.
+        Otherwise it is True.
+        """
+        return any([comp.fit_center for comp in self.components])
+
+    def interpolateStarList(self, stars):
+        """Update the stars to have the current interpolated fit parameters according to the
+        current PSF model.
+
+        :param stars:       List of Star instances to update.
+
+        :returns:           List of Star instances with their fit parameters updated.
+        """
+        stars = setup_params(stars)
+        for k, comp in enumerate(self.components):
+            stars = comp.interpolateStarList(stars)
+        return stars
+
+    def interpolateStar(self, star):
+        """Update the star to have the current interpolated fit parameters according to the
+        current PSF model.
+
+        :param star:        Star instance to update.
+
+        :returns:           Star instance with its fit parameters updated.
+        """
+        star, = self.setup_params([star])
+        for k, comp in enumerate(self.components):
+            star = comp.interpolateStar(star)
+        return star
+
+    def _drawStar(self, star, center=None):
+        # Draw each component
+        comp_stars = []
+        for k, comp in enumerate(self.components):
+            comp_star = comp._drawStar(star, center=center)
+            comp_stars.append(comp_star)
+
+        # Add them up.
+        image = star.image.copy()
+        image.setZero()
+        for comp_star in comp_stars:
+            image += comp_star.image
+        return Star(star.data.withNew(image=image), star.fit)
+
+    def _getRawProfile(self, star):
+        # Get each component profile
+        profiles = []
+        methods = []
+        for k, comp in enumerate(self.components):
+            prof, method = comp._getRawProfile(star)
+            profiles.append(prof)
+            methods.append(method)
+        # For sum, all methods must be the same.  TODO: enforce this a creation.
+        assert all([m == methods[0] for m in methods])
+
+        # Add them up.
+        return galsim.Sum(profiles), methods[0]
+
+    def _finish_write(self, fits, extname, logger):
+        """Finish the writing process with any class-specific steps.
+
+        :param fits:        An open fitsio.FITS object
+        :param extname:     The base name of the extension to write to.
+        :param logger:      A logger object for logging debug info.
+        """
+        logger = galsim.config.LoggerWrapper(logger)
+        chisq_dict = {
+            'chisq' : self.chisq,
+            'last_delta_chisq' : self.last_delta_chisq,
+            'dof' : self.dof,
+            'nremoved' : self.nremoved,
+        }
+        write_kwargs(fits, extname + '_chisq', chisq_dict)
+        logger.debug("Wrote the chisq info to extension %s",extname + '_chisq')
+        for k, comp in enumerate(self.components):
+            comp._write(fits, extname + '_' + str(k), logger=logger)
+        if self.outliers:
+            self.outliers.write(fits, extname + '_outliers')
+            logger.debug("Wrote the PSF outliers to extension %s",extname + '_outliers')
+
+    def _finish_read(self, fits, extname, logger):
+        """Finish the reading process with any class-specific steps.
+
+        :param fits:        An open fitsio.FITS object
+        :param extname:     The base name of the extension to write to.
+        :param logger:      A logger object for logging debug info.
+        """
+        chisq_dict = read_kwargs(fits, extname + '_chisq')
+        for key in chisq_dict:
+            setattr(self, key, chisq_dict[key])
+
+        ncomponents = self.components
+        self.components = []
+        for k in range(ncomponents):
+            self.components.append(PSF._read(fits, extname + '_' + str(k), logger=logger))
+        if extname + '_outliers' in fits:
+            self.outliers = Outliers.read(fits, extname + '_outliers')
+        else:
+            self.outliers = None
+        # Set up all the num's properly now that everything is constructed.
+        self.set_num(None)
diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
new file mode 100644
index 00000000..fe5c7ca2
--- /dev/null
+++ b/tests/test_sumpsf.py
@@ -0,0 +1,147 @@
+# Copyright (c) 2016 by Mike Jarvis and the other collaborators on GitHub at
+# https://github.com/rmjarvis/Piff  All rights reserved.
+#
+# Piff is free software: Redistribution and use in source and binary forms
+# with or without modification, are permitted provided that the following
+# conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the disclaimer given in the accompanying LICENSE
+#    file.
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the disclaimer given in the documentation
+#    and/or other materials provided with the distribution.
+
+import galsim
+import numpy as np
+import piff
+import os
+import fitsio
+
+from piff_test_helper import timer
+
+@timer
+def test_trivial_sum():
+    """Test the trivial case of using a Sum of 1 component.
+    """
+    # This is essentially the same as test_single_image in test_simple.py
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_trivial_sum.log')
+
+    # Make the image
+    image = galsim.Image(2048, 2048, scale=0.26)
+
+    # Where to put the stars.
+    x_list = [ 123.12, 345.98, 567.25, 1094.94, 924.15, 1532.74, 1743.11, 888.39, 1033.29, 1409.31 ]
+    y_list = [ 345.43, 567.45, 1094.32, 924.29, 1532.92, 1743.83, 888.83, 1033.19, 1409.20, 123.11 ]
+
+    # Draw a Gaussian PSF at each location on the image.
+    sigma = 1.3
+    g1 = 0.23
+    g2 = -0.17
+    psf = galsim.Gaussian(sigma=sigma).shear(g1=g1, g2=g2)
+    targets = []
+    for x,y in zip(x_list, y_list):
+        bounds = galsim.BoundsI(int(x-31), int(x+32), int(y-31), int(y+32))
+        psf.drawImage(image=image[bounds], method='no_pixel', center=(x,y))
+        targets.append(image[bounds])
+    image.addNoise(galsim.GaussianNoise(rng=galsim.BaseDeviate(1234), sigma=1e-6))
+
+    # Write out the image to a file
+    image_file = os.path.join('output','trivial_sum_im.fits')
+    image.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x_list), dtype=dtype)
+    data['x'] = x_list
+    data['y'] = y_list
+    cat_file = os.path.join('output','trivial_sum_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    psf_file = os.path.join('output','trivial_sum_psf.fits')
+    stamp_size = 48
+    config = {
+        'input' : {
+            'image_file_name' : image_file,
+            'cat_file_name' : cat_file,
+            'stamp_size' : stamp_size
+        },
+        'psf' : {
+            'type' : 'Sum',
+            'components': [
+                {
+                    'type' : 'Simple',
+                    'model' : { 'type' : 'Gaussian',
+                                'fastfit': True,
+                                'include_pixel': False },
+                    'interp' : { 'type' : 'Mean' },
+                }
+            ],
+            'max_iter' : 10,
+            'chisq_thresh' : 0.2,
+        },
+        'output' : { 'file_name' : psf_file },
+    }
+
+    piff.piffify(config, logger)
+    psf = piff.read(psf_file)
+
+    assert type(psf) is piff.SumPSF
+    assert len(psf.components) == 1
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.Gaussian
+    assert type(psf.components[0].interp) is piff.Mean
+
+    assert psf.chisq_thresh == 0.2
+    assert psf.max_iter == 10
+
+    for i, star in enumerate(psf.stars):
+        target = targets[i]
+        test_star = psf.drawStar(star)
+        true_params = [ sigma, g1, g2 ]
+        np.testing.assert_almost_equal(test_star.fit.params[0], true_params, decimal=4)
+        # The drawn image should be reasonably close to the target.
+        target = target[star.image.bounds]
+        print('target max diff = ',np.max(np.abs(test_star.image.array - target.array)))
+        np.testing.assert_allclose(test_star.image.array, target.array, atol=1.e-5)
+
+        # test that draw works
+        test_image = psf.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                              flux=star.fit.flux, offset=star.fit.center/image.scale)
+        # This image should be very close to the same values as test_star
+        # Just make sure to compare over the same bounds.
+        b = test_star.image.bounds & test_image.bounds
+        print('draw/drawStar max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
+        np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-9)
+
+    # Outlier is a valid option
+    config['psf']['outliers'] = {
+        'type' : 'Chisq',
+        'nsigma' : 5,
+        'max_remove' : 3,
+    }
+    piff.piffify(config)
+    psf2 = piff.read(psf_file)
+    # This didn't remove anything, so the result is the same.
+    # (Use the last star from the above loop for comparison.)
+    test_star2 = psf2.drawStar(star)
+    np.testing.assert_almost_equal(test_star2.fit.params[0], true_params, decimal=4)
+    assert test_star2.image == test_star.image
+
+    test_image2 = psf2.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                           flux=star.fit.flux, offset=star.fit.center/image.scale)
+    assert test_image2 == test_image
+
+    # Error if components is missing
+    config1 = config.copy()
+    del config1['psf']['components']
+    with np.testing.assert_raises(ValueError):
+        piff.process(config)
+
+
+if __name__ == '__main__':
+    test_trivial_sum()

From 9d0e7e9c9e86d5aebc9457bec019d5acfb1b2a61 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Thu, 11 May 2023 10:32:49 -0400
Subject: [PATCH 02/12] Add test of two components

---
 piff/gsobject_model.py |   2 +-
 piff/sumpsf.py         |   1 +
 tests/test_sumpsf.py   | 146 +++++++++++++++++++++++++++++++++++++++--
 3 files changed, 141 insertions(+), 8 deletions(-)

diff --git a/piff/gsobject_model.py b/piff/gsobject_model.py
index 17d23584..78ef7022 100644
--- a/piff/gsobject_model.py
+++ b/piff/gsobject_model.py
@@ -314,7 +314,7 @@ def fit(self, star, fastfit=None, logger=None, convert_func=None):
         if self._centered:
             center = (du, dv)
         else:
-            center = (0.0, 0.0)
+            center = star.fit.center
             params = [du, dv] + params
             params_var = np.concatenate([var[1:3], params_var])
         if self._fit_flux:
diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index 2bd9f4a5..c09fb4ee 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -178,6 +178,7 @@ def single_iteration(self, stars, logger, convert_func):
 
         # Fit each component in order
         for k, comp in enumerate(self.components):
+            logger.info("Starting work on component %d (%s)", k, comp._type_name)
             # Update stars to subtract current fit from other components.
             modified_stars = []
             for i, star in enumerate(stars):
diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
index fe5c7ca2..e3faa5a4 100644
--- a/tests/test_sumpsf.py
+++ b/tests/test_sumpsf.py
@@ -21,7 +21,7 @@
 from piff_test_helper import timer
 
 @timer
-def test_trivial_sum():
+def test_trivial_sum1():
     """Test the trivial case of using a Sum of 1 component.
     """
     # This is essentially the same as test_single_image in test_simple.py
@@ -29,12 +29,12 @@ def test_trivial_sum():
     if __name__ == '__main__':
         logger = piff.config.setup_logger(verbose=2)
     else:
-        logger = piff.config.setup_logger(log_file='output/test_trivial_sum.log')
+        logger = piff.config.setup_logger(log_file='output/test_trivial_sum1.log')
 
     # Make the image
     image = galsim.Image(2048, 2048, scale=0.26)
 
-    # Where to put the stars.
+    # Where to put the stars. (Randomish, but make sure no blending.)
     x_list = [ 123.12, 345.98, 567.25, 1094.94, 924.15, 1532.74, 1743.11, 888.39, 1033.29, 1409.31 ]
     y_list = [ 345.43, 567.45, 1094.32, 924.29, 1532.92, 1743.83, 888.83, 1033.19, 1409.20, 123.11 ]
 
@@ -51,7 +51,7 @@ def test_trivial_sum():
     image.addNoise(galsim.GaussianNoise(rng=galsim.BaseDeviate(1234), sigma=1e-6))
 
     # Write out the image to a file
-    image_file = os.path.join('output','trivial_sum_im.fits')
+    image_file = os.path.join('output','trivial_sum1_im.fits')
     image.write(image_file)
 
     # Write out the catalog to a file
@@ -59,10 +59,10 @@ def test_trivial_sum():
     data = np.empty(len(x_list), dtype=dtype)
     data['x'] = x_list
     data['y'] = y_list
-    cat_file = os.path.join('output','trivial_sum_cat.fits')
+    cat_file = os.path.join('output','trivial_sum1_cat.fits')
     fitsio.write(cat_file, data, clobber=True)
 
-    psf_file = os.path.join('output','trivial_sum_psf.fits')
+    psf_file = os.path.join('output','trivial_sum1_psf.fits')
     stamp_size = 48
     config = {
         'input' : {
@@ -142,6 +142,138 @@ def test_trivial_sum():
     with np.testing.assert_raises(ValueError):
         piff.process(config)
 
+@timer
+def test_easy_sum2():
+    """Test a fairly easy case with 2 components, one constant, the other linear.
+    """
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_eaay_sum2.log')
+
+    # Make the image
+    image = galsim.Image(2048, 2048, scale=0.26)
+
+    # Where to put the stars. (Randomish, but make sure no blending.)
+    x_list = [ 123.12, 345.98, 567.25, 1094.94, 924.15, 1532.74, 1743.11, 888.39, 1033.29, 1409.31 ]
+    y_list = [ 345.43, 567.45, 1094.32, 924.29, 1532.92, 1743.83, 888.83, 1033.19, 1409.20, 123.11 ]
+
+    # The model is a sum of two Gaussians with very different size and interpolation patterns.
+    # A wide Gaussian is the same everywhere
+    # A narrow Gaussian has a linear slope in its parameters across the field of view.
+    sigma1 = 2.3
+    g1 = 0.23
+    g2 = -0.17
+    psf1 = galsim.Gaussian(sigma=sigma1).shear(g1=g1, g2=g2)
+    targets = []
+    for x,y in zip(x_list, y_list):
+        # sigma ranges from 1.2 to 1.5
+        sigma2 = 1.2 + 0.3 * x / 2048
+        # dy ranges from -0.7 arcsec to +0.8 arcsec
+        dy = -0.7 + 1.5 * y / 2048
+        psf2 = galsim.Gaussian(sigma=sigma2).shift(0,dy) * 0.2
+        bounds = galsim.BoundsI(int(x-33), int(x+33), int(y-33), int(y+33))
+        psf = psf1 + psf2
+        psf.drawImage(image=image[bounds], center=(x,y))
+        targets.append(image[bounds])
+    image.addNoise(galsim.GaussianNoise(rng=galsim.BaseDeviate(1234), sigma=1e-6))
+
+    # Write out the image to a file
+    image_file = os.path.join('output','easy_sum2_im.fits')
+    image.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x_list), dtype=dtype)
+    data['x'] = x_list
+    data['y'] = y_list
+    cat_file = os.path.join('output','easy_sum2_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    psf_file = os.path.join('output','easy_sum2_psf.fits')
+    stamp_size = 48
+    config = {
+        'input' : {
+            'image_file_name' : image_file,
+            'cat_file_name' : cat_file,
+            'stamp_size' : stamp_size,
+            'noise': 1.e-12,
+            'trust_pos': True,
+        },
+        'select' : {
+            'max_snr': 1.e6,
+        },
+        'psf' : {
+            'type' : 'Sum',
+            'components': [
+                {
+                    'model' : { 'type' : 'Gaussian',
+                                'include_pixel': True,
+                                'fastfit': True,
+                              },
+                    'interp' : { 'type' : 'Mean' },
+                },
+                {
+                    'model' : { 'type' : 'Gaussian',
+                                'include_pixel': True,
+                                'init': 'zero',
+                                'centered': False,
+                                'fit_flux': True,
+                                'fastfit': True,
+                              },
+                    'interp' : { 'type' : 'Polynomial',
+                                 'order': 1 },
+                }
+            ],
+        },
+        'output' : { 'file_name' : psf_file },
+    }
+
+    piff.piffify(config, logger)
+    psf = piff.read(psf_file)
+
+    assert type(psf) is piff.SumPSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.Gaussian
+    assert type(psf.components[0].interp) is piff.Mean
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.Gaussian
+    assert type(psf.components[1].interp) is piff.Polynomial
+
+    # Defaults
+    assert psf.chisq_thresh == 0.1
+    assert psf.max_iter == 30
+
+    for i, star in enumerate(psf.stars):
+        print('star',i)
+        target = targets[i]
+        test_star = psf.drawStar(star)
+        true_params1 = [ sigma1, g1, g2 ]
+        print('comp0: ',test_star.fit.params[0],' =? ',true_params1)
+        np.testing.assert_almost_equal(test_star.fit.params[0], true_params1, decimal=2)
+        true_params2 = [ 0.2, 0, -0.7 + 1.5*y_list[i]/2048, 1.2 + 0.3*x_list[i]/2048, 0, 0 ]
+        print('comp1: ',test_star.fit.params[1],' =? ',true_params2)
+        np.testing.assert_almost_equal(test_star.fit.params[1], true_params2, decimal=2)
+
+        # The drawn image should be reasonably close to the target.
+        target = target[star.image.bounds]
+        print('target max diff = ',np.max(np.abs(test_star.image.array - target.array)))
+        target.write('junk1.fits')
+        test_star.image.write('junk2.fits')
+        (target-test_star.image).write('junk3.fits')
+        np.testing.assert_allclose(test_star.image.array, target.array, atol=1.e-5)
+
+        # Draw should produce something almost identical (modulo bounds).
+        test_image = psf.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                              flux=star.fit.flux, offset=star.fit.center/image.scale)
+        b = test_star.image.bounds & test_image.bounds
+        print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
+        np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-8)
+
+
 
 if __name__ == '__main__':
-    test_trivial_sum()
+    test_trivial_sum1()
+    test_easy_sum2()

From fe17e308e68ecae7a82f0a5f78e80fad9ebbbd79 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Thu, 11 May 2023 14:48:20 -0400
Subject: [PATCH 03/12] Add test with mixed auto and no_pixel components

---
 piff/sumpsf.py       |  19 ++++++-
 tests/test_sumpsf.py | 130 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 146 insertions(+), 3 deletions(-)

diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index c09fb4ee..f0558f9f 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -262,11 +262,24 @@ def _getRawProfile(self, star):
             prof, method = comp._getRawProfile(star)
             profiles.append(prof)
             methods.append(method)
-        # For sum, all methods must be the same.  TODO: enforce this a creation.
-        assert all([m == methods[0] for m in methods])
+        if any([m != methods[0] for m in methods]):
+            assert all([m == 'no_pixel' or m == 'auto' for m in methods])
+            # Then some of these are auto, and others are no_pixel.
+            # Call the whole thing no_pixel, and manually convolve the auto components.
+            pixel = star.data.local_wcs.toWorld(galsim.Pixel(1.0))
+            new_profiles = []
+            for p,m in zip(profiles, methods):
+                if m == 'auto':
+                    new_profiles.append(galsim.Convolve(p, pixel))
+                else:
+                    new_profiles.append(p)
+            profiles = new_profiles
+            method = 'no_pixel'
+        else:
+            method = methods[0]
 
         # Add them up.
-        return galsim.Sum(profiles), methods[0]
+        return galsim.Sum(profiles), method
 
     def _finish_write(self, fits, extname, logger):
         """Finish the writing process with any class-specific steps.
diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
index e3faa5a4..c2d60824 100644
--- a/tests/test_sumpsf.py
+++ b/tests/test_sumpsf.py
@@ -273,7 +273,137 @@ def test_easy_sum2():
         np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-8)
 
 
+@timer
+def test_mixed_pixel_sum2():
+    """Test case one component including pixel, but not other one.
+    """
+    # This is identical to test_trivial_sum2, except that component 2 is drawn and fit with
+    # no pixel convolution.
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_mixel_pixel.log')
+        logger = piff.config.setup_logger(verbose=2)
+
+    # Make the image
+    image = galsim.Image(2048, 2048, scale=0.26)
+
+    # Where to put the stars.
+    x_list = [ 123.12, 345.98, 567.25, 1094.94, 924.15, 1532.74, 1743.11, 888.39, 1033.29, 1409.31 ]
+    y_list = [ 345.43, 567.45, 1094.32, 924.29, 1532.92, 1743.83, 888.83, 1033.19, 1409.20, 123.11 ]
+
+    # The model is a sum of two Gaussians with very different size and interpolation patterns.
+    # A wide Gaussian is the same everywhere
+    # A narrow Gaussian has a linear slope in its parameters across the field of view.
+    sigma1 = 2.3
+    g1 = 0.23
+    g2 = -0.17
+    psf1 = galsim.Gaussian(sigma=sigma1).shear(g1=g1, g2=g2)
+    targets = []
+    for x,y in zip(x_list, y_list):
+        # sigma ranges from 0.5 to 0.8
+        sigma2 = 0.5 + 0.3 * x / 2048
+        # dy ranges from -0.3 arcsec to +0.2 arcsec
+        dy = -0.3 + 0.5 * y / 2048
+        psf2 = galsim.Gaussian(sigma=sigma2).shift(0,dy) * 0.13
+        bounds = galsim.BoundsI(int(x-33), int(x+33), int(y-33), int(y+33))
+        psf1.drawImage(image=image[bounds], center=(x,y))
+        psf2.drawImage(image=image[bounds], center=(x,y), add_to_image=True, method='no_pixel')
+        targets.append(image[bounds])
+    image.addNoise(galsim.GaussianNoise(rng=galsim.BaseDeviate(1234), sigma=1e-6))
+
+    # Write out the image to a file
+    image_file = os.path.join('output','mixed_pixel_im.fits')
+    image.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x_list), dtype=dtype)
+    data['x'] = x_list
+    data['y'] = y_list
+    cat_file = os.path.join('output','mixed_pixel_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    psf_file = os.path.join('output','mixed_pixel_psf.fits')
+    stamp_size = 48
+    config = {
+        'input' : {
+            'image_file_name' : image_file,
+            'cat_file_name' : cat_file,
+            'stamp_size' : stamp_size,
+            'noise': 1.e-12,
+            'trust_pos': True,
+        },
+        'select' : {
+            'max_snr': 1.e6,
+        },
+        'psf' : {
+            'type' : 'Sum',
+            'components': [
+                {
+                    'model' : { 'type' : 'Gaussian',
+                                'include_pixel': True,
+                                'fastfit': True,
+                              },
+                    'interp' : { 'type' : 'Mean' },
+                },
+                {
+                    'model' : { 'type' : 'Gaussian',
+                                'include_pixel': False,
+                                'init': '(0.1, 0.2)',
+                                'centered': False,
+                                'fit_flux': True,
+                                'fastfit': True,
+                              },
+                    'interp' : { 'type' : 'Polynomial',
+                                 'order': 1 },
+                }
+            ],
+        },
+        'output' : { 'file_name' : psf_file },
+    }
+
+    piff.piffify(config, logger)
+    psf = piff.read(psf_file)
+
+    assert type(psf) is piff.SumPSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.Gaussian
+    assert type(psf.components[0].interp) is piff.Mean
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.Gaussian
+    assert type(psf.components[1].interp) is piff.Polynomial
+
+    for i, star in enumerate(psf.stars):
+        print('star',i)
+        target = targets[i]
+        test_star = psf.drawStar(star)
+        true_params1 = [ sigma1, g1, g2 ]
+        print('comp0: ',test_star.fit.params[0],' =? ',true_params1)
+        np.testing.assert_almost_equal(test_star.fit.params[0], true_params1, decimal=2)
+        true_params2 = [ 0.13, 0, -0.3 + 0.5*y_list[i]/2048, 0.5 + 0.3*x_list[i]/2048, 0, 0 ]
+        print('comp1: ',test_star.fit.params[1],' =? ',true_params2)
+        np.testing.assert_almost_equal(test_star.fit.params[1], true_params2, decimal=2)
+
+        # The drawn image should be reasonably close to the target.
+        target = target[star.image.bounds]
+        print('target max diff = ',np.max(np.abs(test_star.image.array - target.array)))
+        target.write('junk1.fits')
+        test_star.image.write('junk2.fits')
+        (target-test_star.image).write('junk3.fits')
+        np.testing.assert_allclose(test_star.image.array, target.array, atol=1.e-5)
+
+        # Draw should produce something almost identical (modulo bounds).
+        test_image = psf.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                              flux=star.fit.flux, offset=star.fit.center/image.scale)
+        b = test_star.image.bounds & test_image.bounds
+        print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
+        np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-7)
+
 
 if __name__ == '__main__':
     test_trivial_sum1()
     test_easy_sum2()
+    test_mixed_pixel_sum2()

From 893887fcd583ae8f1ef02b5c166f3dc86419bf55 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Tue, 16 May 2023 14:22:11 -0400
Subject: [PATCH 04/12] Add tests with PixelGrid

---
 tests/test_sumpsf.py | 108 ++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 102 insertions(+), 6 deletions(-)

diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
index c2d60824..3163eb06 100644
--- a/tests/test_sumpsf.py
+++ b/tests/test_sumpsf.py
@@ -260,9 +260,6 @@ def test_easy_sum2():
         # The drawn image should be reasonably close to the target.
         target = target[star.image.bounds]
         print('target max diff = ',np.max(np.abs(test_star.image.array - target.array)))
-        target.write('junk1.fits')
-        test_star.image.write('junk2.fits')
-        (target-test_star.image).write('junk3.fits')
         np.testing.assert_allclose(test_star.image.array, target.array, atol=1.e-5)
 
         # Draw should produce something almost identical (modulo bounds).
@@ -272,6 +269,65 @@ def test_easy_sum2():
         print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
         np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-8)
 
+    # Repeat with both components being PixelGrid mdoels.
+    # Both pretty chunky, so this test doesn't take forever.
+    config['psf']['components'][0]['model'] = {
+        'type': 'PixelGrid',
+        'scale': 1.04,  # 4x4 native pixels per grid pixel.
+        'size': 10,     # Covers 40x40 original pixels.
+    }
+    if __name__ == '__main__':
+        # The small Gaussian moves around the central ~8 pixels.
+        # Make the grid big enough to capture the whole moving Gaussian.
+        grid_size = 14
+        tol = 3.e-4
+    else:
+        # For faster running on CI, use a bit smaller central piece with higher tolerance.
+        grid_size = 10
+        tol = 5.e-4
+    config['psf']['components'][1]['model'] = {
+        'type': 'PixelGrid',
+        'scale': 0.26,  # native pixel scale.
+        'size': grid_size,
+        'init': 'zero',
+        'centered': False,
+        'fit_flux': True,
+    }
+    psf = piff.process(config, logger)
+    assert type(psf) is piff.SumPSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.PixelGrid
+    assert type(psf.components[0].interp) is piff.Mean
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.PixelGrid
+    assert type(psf.components[1].interp) is piff.Polynomial
+
+    for i, star in enumerate(psf.stars):
+        print('star',i)
+        target = targets[i]
+        test_star = psf.drawStar(star)
+        # Note: There is no useful test to make of the component parameters.
+        #print('comp0: ',test_star.fit.params[0])
+        #print('comp1: ',test_star.fit.params[1])
+
+        # The drawn image should be reasonably close to the target.
+        # Although not nearly as close as above.
+        # Also, shrink the bounds a bit, since the PixelGrid edge treatment isn't great.
+        b = star.image.bounds.withBorder(-8)
+        target = target[b]
+        print('target max diff = ',np.max(np.abs(test_star.image[b].array - target[b].array)))
+        print('target max value = ',np.max(np.abs(target[b].array)))
+        np.testing.assert_allclose(test_star.image[b].array, target[b].array, atol=tol)
+
+        # Draw should produce something almost identical (modulo bounds).
+        test_image = psf.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                              flux=star.fit.flux, offset=np.array(star.fit.center)/image.scale)
+        b = test_star.image.bounds & test_image.bounds
+        print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
+        np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-7)
+
+
 
 @timer
 def test_mixed_pixel_sum2():
@@ -390,9 +446,6 @@ def test_mixed_pixel_sum2():
         # The drawn image should be reasonably close to the target.
         target = target[star.image.bounds]
         print('target max diff = ',np.max(np.abs(test_star.image.array - target.array)))
-        target.write('junk1.fits')
-        test_star.image.write('junk2.fits')
-        (target-test_star.image).write('junk3.fits')
         np.testing.assert_allclose(test_star.image.array, target.array, atol=1.e-5)
 
         # Draw should produce something almost identical (modulo bounds).
@@ -403,6 +456,49 @@ def test_mixed_pixel_sum2():
         np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-7)
 
 
+    # Repeat with the first component being a PixelGrid.
+    config['psf']['components'][0]['model'] = {
+        'type': 'PixelGrid',
+        'scale': 1.04,  # 4x4 pixel grid
+        'size': 10,
+    }
+    psf = piff.process(config, logger)
+    assert type(psf) is piff.SumPSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.PixelGrid
+    assert type(psf.components[0].interp) is piff.Mean
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.Gaussian
+    assert type(psf.components[1].interp) is piff.Polynomial
+
+    for i, star in enumerate(psf.stars):
+        print('star',i)
+        target = targets[i]
+        test_star = psf.drawStar(star)
+        #print('comp0: ',test_star.fit.params[0])
+        true_params2 = [ 0.13, 0, -0.3 + 0.5*y_list[i]/2048, 0.5 + 0.3*x_list[i]/2048, 0, 0 ]
+        print('comp1: ',test_star.fit.params[1],' =? ',true_params2)
+        # Note: don't test these for exact match anymore.  The only thing we really care about
+        # here is that the drawn image is reasonably close.
+
+        # The drawn image should be reasonably close to the target.
+        # Although not nearly as close as above.
+        # Also, shrink the bounds a bit, since the PixelGrid edge treatment isn't great.
+        b = star.image.bounds.withBorder(-8)
+        target = target[b]
+        print('target max diff = ',np.max(np.abs(test_star.image[b].array - target[b].array)))
+        print('target max value = ',np.max(np.abs(target[b].array)))
+        np.testing.assert_allclose(test_star.image[b].array, target[b].array, atol=3.e-4)
+
+        # Draw should produce something almost identical (modulo bounds).
+        test_image = psf.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                              flux=star.fit.flux, offset=np.array(star.fit.center)/image.scale)
+        b = test_star.image.bounds & test_image.bounds
+        print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
+        np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-7)
+
+
 if __name__ == '__main__':
     test_trivial_sum1()
     test_easy_sum2()

From 4bdf8a04360ec710f3000000f667a6ec65866722 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Thu, 6 Jul 2023 21:11:25 -0400
Subject: [PATCH 05/12] typo

---
 tests/test_sumpsf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
index 3163eb06..1e954b0a 100644
--- a/tests/test_sumpsf.py
+++ b/tests/test_sumpsf.py
@@ -333,7 +333,7 @@ def test_easy_sum2():
 def test_mixed_pixel_sum2():
     """Test case one component including pixel, but not other one.
     """
-    # This is identical to test_trivial_sum2, except that component 2 is drawn and fit with
+    # This is identical to test_easy_sum2, except that component 2 is drawn and fit with
     # no pixel convolution.
 
     if __name__ == '__main__':

From 86930682018057ea663f418f9aaa5b790933cd22 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Thu, 6 Jul 2023 22:47:41 -0400
Subject: [PATCH 06/12] coverage

---
 piff/sumpsf.py       |  2 +-
 tests/test_sumpsf.py | 12 +++++++-----
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index f0558f9f..11821012 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -222,7 +222,7 @@ def interpolateStarList(self, stars):
 
         :returns:           List of Star instances with their fit parameters updated.
         """
-        stars = setup_params(stars)
+        stars = self.setup_params(stars)
         for k, comp in enumerate(self.components):
             stars = comp.interpolateStarList(stars)
         return stars
diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
index 1e954b0a..46fc4508 100644
--- a/tests/test_sumpsf.py
+++ b/tests/test_sumpsf.py
@@ -246,10 +246,11 @@ def test_easy_sum2():
     assert psf.chisq_thresh == 0.1
     assert psf.max_iter == 30
 
-    for i, star in enumerate(psf.stars):
+    test_stars = psf.interpolateStarList(psf.stars)
+    test_stars = psf.drawStarList(test_stars)
+    for i, test_star in enumerate(test_stars):
         print('star',i)
         target = targets[i]
-        test_star = psf.drawStar(star)
         true_params1 = [ sigma1, g1, g2 ]
         print('comp0: ',test_star.fit.params[0],' =? ',true_params1)
         np.testing.assert_almost_equal(test_star.fit.params[0], true_params1, decimal=2)
@@ -258,13 +259,14 @@ def test_easy_sum2():
         np.testing.assert_almost_equal(test_star.fit.params[1], true_params2, decimal=2)
 
         # The drawn image should be reasonably close to the target.
-        target = target[star.image.bounds]
+        target = target[test_star.image.bounds]
         print('target max diff = ',np.max(np.abs(test_star.image.array - target.array)))
         np.testing.assert_allclose(test_star.image.array, target.array, atol=1.e-5)
 
         # Draw should produce something almost identical (modulo bounds).
-        test_image = psf.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
-                              flux=star.fit.flux, offset=star.fit.center/image.scale)
+        test_image = psf.draw(x=test_star['x'], y=test_star['y'],
+                              stamp_size=config['input']['stamp_size'],
+                              flux=test_star.fit.flux, offset=test_star.fit.center/image.scale)
         b = test_star.image.bounds & test_image.bounds
         print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
         np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-8)

From 89416479d770767e3540771460b117f27a9942f0 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Tue, 18 Jul 2023 15:33:23 -0400
Subject: [PATCH 07/12] Edits based on Theo's code review.

---
 piff/sumpsf.py       | 16 +++++++---------
 tests/test_sumpsf.py |  2 +-
 2 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index 11821012..fd5e0b88 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -54,7 +54,7 @@ def __init__(self, components, outliers=None, chisq_thresh=0.1, max_iter=30):
         self.chisq_thresh = chisq_thresh
         self.max_iter = max_iter
         self.kwargs = {
-            # If components is a list, mark the number of components here for I/O purposed.
+            # If components is a list, mark the number of components here for I/O purposes.
             # But if it's just a number, leave it alone.
             'components': len(components) if isinstance(components, list) else components,
             'outliers': 0,
@@ -88,10 +88,9 @@ def set_num(self, num):
                 self._num_components += comp.num_components
             self._num = self._nums[0]
         else:
+            # else components are not yet built. This will be called again when that's done.
             self._num = None
 
-        # else components are not yet built. This will be called again when that's done.
-
     @property
     def num_components(self):
         return self._num_components
@@ -172,9 +171,7 @@ def single_iteration(self, stars, logger, convert_func):
         nremoved = 0  # For this iteration
 
         # Draw the current models for each component
-        current_models = []
-        for k, comp in enumerate(self.components):
-            current_models.append(comp.drawStarList(stars))
+        current_models = [comp.drawStarList(stars) for comp in self.components]
 
         # Fit each component in order
         for k, comp in enumerate(self.components):
@@ -192,9 +189,10 @@ def single_iteration(self, stars, logger, convert_func):
             # Fit this component
             new_stars, nremoved1 = comp.single_iteration(modified_stars, logger, convert_func)
 
-            # The fit components in new_stars are what we want, but the data parts are
-            # corrupted by subtracting the other components, which we don't want.
-            # Just copy over the fits.
+            # new_stars now has the updated fit components, but the data parts are corrupted by
+            # subtracting the other components during fitting.
+            # Update the stars by copying the undisturbed data from the original stars and the
+            # updated fits from new_stars.
             stars = [Star(s1.data, s2.fit) for s1,s2 in zip(stars, new_stars)]
             nremoved += nremoved1
 
diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
index 46fc4508..8c6af59e 100644
--- a/tests/test_sumpsf.py
+++ b/tests/test_sumpsf.py
@@ -271,7 +271,7 @@ def test_easy_sum2():
         print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
         np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-8)
 
-    # Repeat with both components being PixelGrid mdoels.
+    # Repeat with both components being PixelGrid models.
     # Both pretty chunky, so this test doesn't take forever.
     config['psf']['components'][0]['model'] = {
         'type': 'PixelGrid',

From 7edfe4d6451f21f5c7290b7455be9fbe4ee22723 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Wed, 27 Mar 2024 12:53:20 -0400
Subject: [PATCH 08/12] Remove some more unnecessary enumerates

---
 piff/sumpsf.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index fd5e0b88..2182209c 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -156,7 +156,7 @@ def initialize_params(self, stars, logger, default_init=None):
         stars = self.setup_params(stars)
 
         # Now initialize all the components
-        for k, comp in enumerate(self.components):
+        for comp in self.components:
             stars, nremoved1 = comp.initialize_params(stars, logger, default_init=default_init)
             nremoved += nremoved1
             # After the first one, set default_init to 'zero'
@@ -221,7 +221,7 @@ def interpolateStarList(self, stars):
         :returns:           List of Star instances with their fit parameters updated.
         """
         stars = self.setup_params(stars)
-        for k, comp in enumerate(self.components):
+        for comp in self.components:
             stars = comp.interpolateStarList(stars)
         return stars
 
@@ -234,14 +234,14 @@ def interpolateStar(self, star):
         :returns:           Star instance with its fit parameters updated.
         """
         star, = self.setup_params([star])
-        for k, comp in enumerate(self.components):
+        for comp in self.components:
             star = comp.interpolateStar(star)
         return star
 
     def _drawStar(self, star, center=None):
         # Draw each component
         comp_stars = []
-        for k, comp in enumerate(self.components):
+        for comp in self.components:
             comp_star = comp._drawStar(star, center=center)
             comp_stars.append(comp_star)
 
@@ -256,7 +256,7 @@ def _getRawProfile(self, star):
         # Get each component profile
         profiles = []
         methods = []
-        for k, comp in enumerate(self.components):
+        for comp in self.components:
             prof, method = comp._getRawProfile(star)
             profiles.append(prof)
             methods.append(method)

From 6180dcd41a69ee1d9b6a678764f66c470050c8b2 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Wed, 27 Mar 2024 17:11:36 -0400
Subject: [PATCH 09/12] Don't use pytest 8

---
 .github/workflows/ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 70f8fbc1..aa1ef20d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -89,7 +89,7 @@ jobs:
                 # Extra packages needed for testing
                 # Note: Pin pillow <10 until this bug is fixed:
                 #       https://github.com/python-pillow/Pillow/issues/7259
-                pip install -U nose coverage pytest nbval ipykernel "pillow<10"
+                pip install -U nose coverage "pytest<8" nbval ipykernel "pillow<10"
 
             - name: Install Pixmappy (not on pip)
               run: |

From d99765841ccfd803e13564a9ec1b6ba5c8b8489f Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Fri, 29 Mar 2024 00:40:44 -0400
Subject: [PATCH 10/12] Add option to include or not include reserve and
 flagged stars in StarImages plot

---
 piff/star_stats.py  | 36 +++++++++++++++++++++++++++++-------
 tests/test_stats.py | 30 +++++++++++++++++++++++++++---
 2 files changed, 56 insertions(+), 10 deletions(-)

diff --git a/piff/star_stats.py b/piff/star_stats.py
index d5989081..51ae49d2 100644
--- a/piff/star_stats.py
+++ b/piff/star_stats.py
@@ -41,15 +41,23 @@ class StarStats(Stats):
                             (without replacement). [default: 10]
     :param adjust_stars:    Boolean. If true, when computing, will also fit for best
                             starfit center and flux to match observed star. [default: False]
+    :param include_reserve: Whether to inlude reserve stars. [default: True]
+    :param only_reserve:    Whether to skip plotting non-reserve stars. [default: False]
+    :param include_flaggede: Whether to include plotting flagged stars. [default: False]
     :param file_name:       Name of the file to output to. [default: None]
     :param logger:          A logger object for logging debug info. [default: None]
     """
     _type_name = 'StarImages'
 
-    def __init__(self, nplot=10, adjust_stars=False, file_name=None, logger=None):
+    def __init__(self, nplot=10, adjust_stars=False,
+                 include_reserve=True, only_reserve=False, include_flagged=False,
+                 file_name=None, logger=None):
         self.nplot = nplot
         self.file_name = file_name
         self.adjust_stars = adjust_stars
+        self.include_reserve = include_reserve
+        self.only_reserve = only_reserve
+        self.include_flagged = include_flagged
 
     def compute(self, psf, stars, logger=None):
         """
@@ -58,12 +66,21 @@ def compute(self, psf, stars, logger=None):
         :param logger:      A logger object for logging debug info. [default: None]
         """
         logger = galsim.config.LoggerWrapper(logger)
-        # get the shapes
+        # Determine which stars to plot
+        possible_indices = []
+        if self.include_reserve:
+            possible_indices += [i for i,s in enumerate(stars)
+                                 if s.is_reserve and (self.include_flagged or not s.is_flagged)]
+        if not self.only_reserve:
+            possible_indices += [i for i,s in enumerate(stars)
+                                 if not s.is_reserve and (self.include_flagged or not s.is_flagged)]
+        possible_indices = sorted(possible_indices)
+
         if self.nplot == 0 or self.nplot >= len(stars):
-            # select all stars
-            self.indices = np.arange(len(stars))
+            # select all viable stars
+            self.indices = possible_indices
         else:
-            self.indices = np.random.choice(len(stars), self.nplot, replace=False)
+            self.indices = np.random.choice(possible_indices, self.nplot, replace=False)
 
         logger.info("Making {0} Model Stars".format(len(self.indices)))
         self.stars = []
@@ -110,8 +127,13 @@ def plot(self, logger=None, **kwargs):
             ii = i // 2
             jj = (i % 2) * 3
 
-            axs[ii][jj+0].set_title('Star {0}'.format(index))
-            axs[ii][jj+1].set_title('PSF at (u,v) = \n ({0:+.02e}, {1:+.02e})'.format(u, v))
+            title = f'Star {index}'
+            if star.is_reserve:
+                title = 'Reserve ' + title
+            if star.is_flagged:
+                title = 'Flagged ' + title
+            axs[ii][jj+0].set_title(title)
+            axs[ii][jj+1].set_title(f'PSF at (u,v) = \n ({u:+.02e}, {v:+.02e})')
             axs[ii][jj+2].set_title('Star - PSF')
 
             star_image = star.image
diff --git a/tests/test_stats.py b/tests/test_stats.py
index 2ba33b78..992e7cd2 100644
--- a/tests/test_stats.py
+++ b/tests/test_stats.py
@@ -509,6 +509,10 @@ def test_starstats_config():
             'cat_file_name' : cat_file,
             'stamp_size' : 48
         },
+        'select' : {
+            'reserve_frac': 0.2,
+            'seed': 12345,
+        },
         'psf' : {
             'model' : { 'type' : 'Gaussian',
                         'fastfit': True,
@@ -547,7 +551,8 @@ def test_starstats_config():
     # check default nplot
     psf = piff.read(psf_file)
     starStats = piff.StarStats()
-    orig_stars, wcs, pointing = piff.Input.process(config['input'], logger)
+    orig_stars, wcs, pointing = piff.Input.process(config['input'], logger=logger)
+    orig_stars = piff.Select.process(config['select'], orig_stars, logger=logger)
     with np.testing.assert_raises(RuntimeError):
         starStats.write()  # Cannot write before compute
     starStats.compute(psf, orig_stars)
@@ -580,8 +585,22 @@ def test_starstats_config():
     np.testing.assert_array_equal(starStats.indices, np.arange(len(orig_stars)))
     starStats.plot()  # Make sure this runs without error.
 
+    # With include_reserve=False, only 8 stars
+    print('All stars: n=',len(starStats.stars))  # 10 stars total
+    assert len(starStats.stars) == 10
+    starStats = piff.StarStats(nplot=0, include_reserve=False)
+    starStats.compute(psf, orig_stars)
+    assert len(starStats.stars) == 8
+    starStats.plot()  # Make sure this runs without error.
+
+    # With only_reserve=True, only 2 stars
+    starStats = piff.StarStats(nplot=0, only_reserve=True)
+    starStats.compute(psf, orig_stars)
+    assert len(starStats.stars) == 2
+    starStats.plot()  # Make sure this runs without error.
+
     # rerun with adjust stars and see if it did the right thing
-    # first with starstats == False
+    # first with adjust_stars == False
     starStats = piff.StarStats(nplot=0, adjust_stars=False)
     starStats.compute(psf, orig_stars, logger=logger)
     fluxs_noadjust = np.array([s.fit.flux for s in starStats.stars])
@@ -591,7 +610,7 @@ def test_starstats_config():
     # check that ds are 0
     np.testing.assert_array_equal(ds_noadjust, 0)
 
-    # now with starstats == True
+    # now with adjust_stars == True
     starStats = piff.StarStats(nplot=0, adjust_stars=True)
     starStats.compute(psf, orig_stars, logger=logger)
     fluxs_adjust = np.array([s.fit.flux for s in starStats.stars])
@@ -773,6 +792,11 @@ def test_bad_hsm():
                     'type': 'StarImages',
                     'file_name': star_file,
                 },
+                {
+                    'type': 'StarImages',
+                    'file_name': star_file,
+                    'include_flagged': True,
+                },
                 {
                     'type': 'SizeMag',
                     'file_name': sizemag_file,

From 72427145aa18185108513c06e6bbca089a0f6b02 Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Fri, 29 Mar 2024 15:58:53 -0400
Subject: [PATCH 11/12] Make sure convert xrange to int, since it isn't int for
 some Interpolants.

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

diff --git a/piff/pixelgrid.py b/piff/pixelgrid.py
index f3ee0b43..3d7a525c 100644
--- a/piff/pixelgrid.py
+++ b/piff/pixelgrid.py
@@ -366,7 +366,7 @@ def chisq(self, star, logger=None, convert_func=None):
                 v = v[mask] / self.scale
                 ui,uf = np.divmod(u,1)
                 vi,vf = np.divmod(v,1)
-                xr = self.interp.xrange
+                xr = int(np.ceil(self.interp.xrange))
                 # Note arguments are basis pixel position minus image pixel position.
                 # Hence the minus sign in front of uf.
                 argu = -uf[:,np.newaxis] + np.arange(-xr+1,xr+1)[np.newaxis,:]

From dafc650092ff908999a028f92601eafae6860abf Mon Sep 17 00:00:00 2001
From: Mike Jarvis <michael@jarvis.net>
Date: Fri, 29 Mar 2024 16:58:19 -0400
Subject: [PATCH 12/12] Add test that uses PixelGrid to clean up optical part
 of PSF

---
 tests/test_sumpsf.py | 207 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 207 insertions(+)

diff --git a/tests/test_sumpsf.py b/tests/test_sumpsf.py
index 8c6af59e..64edf949 100644
--- a/tests/test_sumpsf.py
+++ b/tests/test_sumpsf.py
@@ -500,8 +500,215 @@ def test_mixed_pixel_sum2():
         print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
         np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-7)
 
+@timer
+def test_atmpsf():
+    # Let reality be atmospheric + optical PSF with realistic variation.
+    # Use Kolmogorov with GP for primary model and PixelGrid to mop up the residuals.
+    # The residuals are mostly just in the center of the profile, so the PixelGrid doesn't need
+    # to be very large to fix the main effects of the optical component.
+
+    # Some parameters copied from psf_wf_movie.py in GalSim repo.
+    Ellerbroek_alts = [0.0, 2.58, 5.16, 7.73, 12.89, 15.46]  # km
+    Ellerbroek_weights = [0.652, 0.172, 0.055, 0.025, 0.074, 0.022]
+    Ellerbroek_interp = galsim.LookupTable(Ellerbroek_alts, Ellerbroek_weights,
+                                           interpolant='linear')
+
+    # Use given number of uniformly spaced altitudes
+    alts = np.max(Ellerbroek_alts)*np.arange(6)/5.
+    weights = Ellerbroek_interp(alts)  # interpolate the weights
+    weights /= sum(weights)  # and renormalize
+
+    spd = []  # Wind speed in m/s
+    dirn = [] # Wind direction in radians
+    r0_500 = [] # Fried parameter in m at a wavelength of 500 nm.
+    u = galsim.UniformDeviate(1234)
+    for i in range(6):
+        spd.append(u() * 20)
+        dirn.append(u() * 360*galsim.degrees)
+        r0_500.append(0.1 * weights[i]**(-3./5))
+
+    screens = galsim.Atmosphere(r0_500=r0_500, speed=spd, direction=dirn, altitude=alts, rng=u,
+                                screen_size=102.4, screen_scale=0.1)
+
+    # Add in an optical screen
+    screens.append(galsim.OpticalScreen(diam=8, defocus=0.7, astig1=-0.8, astig2=0.7,
+                                        trefoil1=-0.6, trefoil2=0.5,
+                                        coma1=0.5, coma2=0.7, spher=0.8, obscuration=0.4))
+
+    aper = galsim.Aperture(diam=8, obscuration=0.4)
+
+    if __name__ == '__main__':
+        size = 2048
+        nstars = 200
+    else:
+        size = 1024
+        nstars = 50
+
+    pixel_scale = 0.2
+    im = galsim.ImageF(size, size, scale=pixel_scale)
+
+    x = u.np.uniform(25, size-25, size=nstars)
+    y = u.np.uniform(25, size-25, size=nstars)
+
+    for k in range(nstars):
+        flux = 100000
+        theta = ((x[k] - size/2) * pixel_scale * galsim.arcsec,
+                 (y[k] - size/2) * pixel_scale * galsim.arcsec)
+
+        psf = screens.makePSF(lam=500, aper=aper, exptime=100, flux=flux, theta=theta)
+        psf.drawImage(image=im, center=(x[k],y[k]), method='phot', rng=u, add_to_image=True)
+        bounds = galsim.BoundsI(int(x[k]-33), int(x[k]+33), int(y[k]-33), int(y[k]+33))
+ 
+    # Add a little noise
+    noise = 10
+    im.addNoise(galsim.GaussianNoise(rng=u, sigma=noise))
+    image_file = os.path.join('output', 'atmpsf_im.fits')
+    im.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x), dtype=dtype)
+    data['x'] = x
+    data['y'] = y
+    cat_file = os.path.join('output','atmpsf_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    psf_file = os.path.join('output','atmpsf.fits')
+    stamp_size = 48
+    config = {
+        'input' : {
+            'image_file_name' : image_file,
+            'cat_file_name' : cat_file,
+            'stamp_size' : 32,
+            'noise': noise**2,
+        },
+        'select' : {
+            'max_snr': 1.e6,
+            'max_edge_frac': 0.1,
+            'hsm_size_reject': True,
+        },
+        'psf' : {
+            'type' : 'Sum',
+            'components': [
+                {
+                    'model' : { 'type' : 'Kolmogorov',
+                                'fastfit': True,
+                              },
+                    'interp' : { 'type' : 'GP',
+                               },
+                },
+                {
+                    'model' : { 'type' : 'PixelGrid',
+                                'init': 'zero',
+                                'scale': pixel_scale,
+                                'size': 15,
+                                'fit_flux': True,
+                              },
+                    'interp' : { 'type' : 'Mean', },
+                }
+            ],
+            'outliers': {
+                'type' : 'Chisq',
+                'nsigma' : 5,
+                'max_remove' : 3,
+            }
+        },
+        'output' : {
+            'file_name' : psf_file,
+            'stats' : [
+                {
+                    'type': 'StarImages',
+                    'file_name': os.path.join('output','test_atmpsf_stars.png'),
+                    'nplot': 10,
+                    'adjust_stars': True,
+                }
+            ],
+        },
+    }
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=1)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_atmpsf.log')
+        logger = piff.config.setup_logger(verbose=1)
+
+    psf = piff.process(config, logger)
+
+    if __name__ == '__main__':
+        output = piff.Output.process(config['output'], logger=logger)
+        output.write(psf, logger=logger)
+
+    assert type(psf) is piff.SumPSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.Kolmogorov
+    assert type(psf.components[0].interp) is piff.GPInterp
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.PixelGrid
+    assert type(psf.components[1].interp) is piff.Mean
+
+    mean_max_rel_err = 0
+    mean_mse = 0
+    count = 0
+    for i, star in enumerate(psf.stars):
+        if star.is_flagged:
+            continue
+        test_star = psf.drawStar(star)
+
+        b = star.image.bounds.withBorder(-8)
+        max_diff = np.max(np.abs(test_star.image[b].array - star.image[b].array))
+        max_val = np.max(np.abs(star.image[b].array))
+        #print('max_diff / max_val = ',max_diff, max_val, max_diff / max_val)
+        mse = np.sum((test_star.image[b].array - star.image[b].array)**2) / flux**2
+        #print('mse = ',mse)
+        mean_max_rel_err += max_diff/max_val
+        mean_mse += mse
+        count += 1
+
+    mean_max_rel_err /= count
+    mean_mse /= count
+    print('mean maximum relative error = ',mean_max_rel_err)
+    print('mean mean-squared error = ',mean_mse)
+    assert mean_max_rel_err < 0.04
+    assert mean_mse < 1.2e-5
+
+    # Without the PixelGrid clean up, it's quite a bit worse.
+    config['psf'] = config['psf']['components'][0]
+    psf = piff.process(config, logger)
+
+    if __name__ == '__main__':
+        config['output']['stats'][0]['file_name'] = os.path.join('output', 'test_atmpsf_stars2.png')
+        output = piff.Output.process(config['output'], logger=logger)
+        output.write(psf, logger=logger)
+
+    mean_max_rel_err = 0
+    mean_mse = 0
+    count = 0
+    for i, star in enumerate(psf.stars):
+        if star.is_flagged:
+            continue
+        test_star = psf.drawStar(star)
+
+        b = star.image.bounds.withBorder(-8)
+        max_diff = np.max(np.abs(test_star.image[b].array - star.image[b].array))
+        max_val = np.max(np.abs(star.image[b].array))
+        #print('max_diff / max_val = ',max_diff / max_val)
+        mse = np.sum((test_star.image[b].array - star.image[b].array)**2) / flux**2
+        #print('mse = ',mse)
+        mean_max_rel_err += max_diff/max_val
+        mean_mse += mse
+        count += 1
+
+    mean_max_rel_err /= count
+    mean_mse /= count
+    print('mean maximum relative error = ',mean_max_rel_err)
+    print('mean mean-squared error = ',mean_mse)
+    assert mean_max_rel_err > 0.06
+    assert mean_mse > 6.e-5
+
 
 if __name__ == '__main__':
     test_trivial_sum1()
     test_easy_sum2()
     test_mixed_pixel_sum2()
+    test_atmpsf()