Address comments

include/lsst/afw/detection/Psf.h

@@ -291,11 +291,9 @@ class Psf : public afw::table::io::PersistableFacade<Psf>,
     explicit Psf(bool isFixed = false, std::size_t capacity = 100);
 
     /**
-     *
-     * This protected virtual member was formerly private (see message below).  We switched it to
-     * protected so that the Psf trampoline class can call it if a python-implemented Psf-derived
-     * class doesn't provide an override.  In contrast, pure virtual private members may remain
-     * private since they must be overridden.
+     * This virtual member is protected (rather than private) so that python-implemented derived
+     * classes may opt to use the default implementation.  C++ derived classes may override this
+     * method, but should not call it.
      */
     virtual std::shared_ptr<Image> doComputeImage(lsst::geom::Point2D const& position,
                                                   image::Color const& color) const;

include/lsst/afw/detection/python.h

@@ -28,6 +28,8 @@
 #include "lsst/afw/detection/Psf.h"
 #include "lsst/afw/typehandling/python.h"
 
+using lsst::afw::typehandling::StorableHelper;
+
 namespace lsst {
 namespace afw {
 namespace detection {
@@ -36,7 +38,7 @@ namespace detection {
  * "Trampoline" for Psf to let it be used as a base class in Python.
  *
  * Subclasses of Psf that are wrapped in %pybind11 should have a similar
- * helper that subclasses `PyPsf<subclass>`. This helper can be
+ * helper that subclasses `PsfTrampoline<subclass>`. This helper can be
  * skipped if the subclass neither adds any virtual methods nor implements
  * any abstract methods.
  *
@@ -45,13 +47,30 @@ namespace detection {
  * @see [pybind11 documentation](https://pybind11.readthedocs.io/en/stable/advanced/classes.html)
  */
 template <typename Base = Psf>
-class PyPsf : public lsst::afw::typehandling::StorableHelper<Base> {
+class PsfTrampoline : public StorableHelper<Base> {
 public:
-    using lsst::afw::typehandling::StorableHelper<Base>::StorableHelper;  // Inherit ctors
     using Image = typename Base::Image;
 
+    /**
+     * Delegating constructor for wrapped class.
+     *
+     * While we would like to simply inherit base class constructors, when doing so, we cannot
+     * change their access specifiers.  One consequence is that it's not possible to use inheritance
+     * to expose a protected constructor to python.  The alternative, used here, is to create a new
+     * public constructor that delegates to the base class public or protected constructor with the
+     * same signature.
+     *
+     * @tparam Args  Variadic type specification
+     * @param ...args  Arguments to forward to the Base class constructor.
+     */
+    template<typename... Args>
+    PsfTrampoline<Base>(Args... args) : StorableHelper<Base>(args...) {}
+
     std::shared_ptr<Psf> clone() const override {
-        PYBIND11_OVERLOAD_PURE(std::shared_ptr<Psf>, Base, clone,);
+        /* __deepcopy__ takes an optional dict, but PYBIND11_OVERLOAD_* won't
+         * compile unless you give it arguments that work for the C++ method
+         */
+        PYBIND11_OVERLOAD_PURE_NAME(std::shared_ptr<Psf>, Base, "__deepcopy__", clone,);
     }
 
     std::shared_ptr<Psf> resized(int width, int height) const override {
@@ -63,12 +82,12 @@ class PyPsf : public lsst::afw::typehandling::StorableHelper<Base> {
     }
 
     // Private and protected c++ members are overloaded to python using underscores.
-    std::shared_ptr<typename Base::Image> doComputeImage(
+    std::shared_ptr<Image> doComputeImage(
         lsst::geom::Point2D const& position,
         image::Color const& color
     ) const override {
         PYBIND11_OVERLOAD_NAME(
-            std::shared_ptr<typename Base::Image>, Base, "_doComputeImage", doComputeImage, position, color
+            std::shared_ptr<Image>, Base, "_doComputeImage", doComputeImage, position, color
         );
     }
 

include/lsst/afw/typehandling/python.h

@@ -51,6 +51,21 @@ class StorableHelper : public Base {
 public:
     using Base::Base;
 
+    /**
+     * Delegating constructor for wrapped class.
+     *
+     * While we would like to simply inherit base class constructors, when doing so, we cannot
+     * change their access specifiers.  One consequence is that it's not possible to use inheritance
+     * to expose a protected constructor to python.  The alternative, used here, is to create a new
+     * public constructor that delegates to the base class public or protected constructor with the
+     * same signature.
+     *
+     * @tparam Args  Variadic type specification
+     * @param ...args  Arguments to forward to the Base class constructor.
+     */
+    template<typename... Args>
+    StorableHelper<Base>(Args... args) : Base(args...) {}
+
     std::shared_ptr<Storable> cloneStorable() const override {
         /* __deepcopy__ takes an optional dict, but PYBIND11_OVERLOAD_* won't
          * compile unless you give it arguments that work for the C++ method

python/lsst/afw/detection/_psf.cc

@@ -26,6 +26,7 @@
 #include <pybind11/pybind11.h>
 
 #include "lsst/utils/python.h"
+#include "lsst/utils/python/PySharedPtr.h"
 
 #include "lsst/geom/Point.h"
 #include "lsst/afw/image/Color.h"
@@ -37,6 +38,8 @@
 namespace py = pybind11;
 using namespace pybind11::literals;
 
+using lsst::utils::python::PySharedPtr;
+
 namespace lsst {
 namespace afw {
 namespace detection {
@@ -52,10 +55,12 @@ void wrapPsf(utils::python::WrapperCollection& wrappers) {
     wrappers.addSignatureDependency("lsst.afw.fits");
 
     auto clsPsf = wrappers.wrapType(
-            py::class_<Psf, std::shared_ptr<Psf>, typehandling::Storable, PyPsf<>>(wrappers.module, "Psf"),
+            py::class_<Psf, PySharedPtr<Psf>, typehandling::Storable, PsfTrampoline<>>(
+                wrappers.module, "Psf"
+            ),
             [](auto& mod, auto& cls) {
                 table::io::python::addPersistableMethods<Psf>(cls);
-                cls.def(py::init<>());  // Dummy constructor for pure-Python subclasses
+                cls.def(py::init<bool, size_t>(), "isFixed"_a=false, "capacity"_a=100);  // Constructor for pure-Python subclasses
                 cls.def("clone", &Psf::clone);
                 cls.def("resized", &Psf::resized, "width"_a, "height"_a);
                 cls.def("computeImage", &Psf::computeImage, "position"_a = NullPoint,

tests/SConscript

@@ -1,7 +1,7 @@
 # -*- python -*-
 from lsst.sconsUtils import scripts, env, targets
 
-pybind11_test_modules = ['testTableArchivesLib', 'testGenericMapLib']
+pybind11_test_modules = ['testTableArchivesLib', 'testGenericMapLib', 'testPsfTrampolineLib']
 
 scripts.BasicSConscript.pybind11(pybind11_test_modules, addUnderscore=False)
 

tests/testPsfTrampolineLib.cc

@@ -0,0 +1,59 @@
+/*
+ * This file is part of afw.
+ *
+ * Developed for the LSST Data Management System.
+ * This product includes software developed by the LSST Project
+ * (https://www.lsst.org).
+ * See the COPYRIGHT file at the top-level directory of this distribution
+ * for details of code ownership.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+// Testing helper module that let's us conveniently invoke c++ methods that should
+// get overridden in python-derived subclasses of Psf.
+
+#include <pybind11/pybind11.h>
+
+#include "lsst/afw/detection/Psf.h"
+
+namespace lsst {
+namespace afw {
+namespace detection {
+
+namespace {
+
+std::shared_ptr<Psf> resizedPsf(const Psf& psf, int nx, int ny) {
+    return psf.resized(nx, ny);
+}
+
+std::shared_ptr<Psf> clonedPsf(const Psf& psf) {
+    return psf.clone();
+}
+
+std::shared_ptr<typehandling::Storable> clonedStorablePsf(const Psf& psf) {
+    return psf.cloneStorable();
+}
+
+}
+
+PYBIND11_MODULE(testPsfTrampolineLib, mod) {
+    mod.def("resizedPsf", &resizedPsf);
+    mod.def("clonedPsf", &clonedPsf);
+    mod.def("clonedStorablePsf", &clonedStorablePsf);
+}
+
+}  // namespace detection
+}  // namespace afw
+}  // namespace lsst

tests/test_psf_trampoline.py

@@ -20,26 +20,30 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import unittest
+from copy import deepcopy
 
 import numpy as np
 
 import lsst.utils.tests
 from lsst.afw.detection import Psf, GaussianPsf
 from lsst.afw.image import Image
-from lsst.geom import Box2I, Point2I, Extent2I
+from lsst.geom import Box2I, Extent2I, Point2I, Point2D
 from lsst.afw.geom.ellipses import Quadrupole
+import testPsfTrampolineLib as cppLib
 
 
 # Subclass Psf in python.  Main tests here are that python virtual methods get
-# resolved by trampoline class.
+# resolved by trampoline class.  The test suite below calls python compute*
+# methods which are implemented in c++ to call the _doCompute* methods defined
+# in the PyGaussianPsf class.
 class PyGaussianPsf(Psf):
     def __init__(self, width, height, sigma):
-        Psf.__init__(self)
+        Psf.__init__(self, isFixed=True)
         self.dimensions = Extent2I(width, height)
         self.sigma = sigma
 
     # "public" virtual overrides
-    def clone(self):
+    def __deepcopy__(self, memo=None):
         return PyGaussianPsf(self.dimensions.x, self.dimensions.y, self.sigma)
 
     def resized(self, width, height):
@@ -65,68 +69,139 @@ def _doComputeApertureFlux(self, radius, position=None, color=None):
         return 1 - np.exp(-0.5*(radius/self.sigma)**2)
 
 
-class PythonPsfInheritTestSuite(lsst.utils.tests.TestCase):
-    def testTrampoline(self):
-        params = [
+class PsfTrampolineTestSuite(lsst.utils.tests.TestCase):
+    def setUp(self):
+        self.pgps = []
+        self.gps = []
+        for width, height, sigma in [
             (5, 5, 1.1),
             (5, 3, 1.2),
             (7, 7, 1.3)
-        ]
-        for width, height, sigma in params:
-            pgp = PyGaussianPsf(width, height, sigma)
-            gp = GaussianPsf(width, height, sigma)
+        ]:
+            self.pgps.append(PyGaussianPsf(width, height, sigma))
+            self.gps.append(GaussianPsf(width, height, sigma))
 
+    def testImages(self):
+        for pgp, gp in zip(self.pgps, self.gps):
             self.assertImagesAlmostEqual(
                 pgp.computeImage(),
                 gp.computeImage()
             )
-
             self.assertImagesAlmostEqual(
                 pgp.computeKernelImage(),
                 gp.computeKernelImage()
             )
 
+    def testApertureFlux(self):
+        for pgp, gp in zip(self.pgps, self.gps):
             for r in [0.1, 0.2, 0.3]:
                 self.assertAlmostEqual(
                     pgp.computeApertureFlux(r),
                     gp.computeApertureFlux(r)
                 )
 
+    def testPeak(self):
+        for pgp, gp in zip(self.pgps, self.gps):
             self.assertAlmostEqual(
                 pgp.computePeak(),
                 gp.computePeak()
             )
 
+    def testBBox(self):
+        for pgp, gp in zip(self.pgps, self.gps):
             self.assertEqual(
                 pgp.computeBBox(),
                 gp.computeBBox()
             )
 
-            self.assertEqual(
+    def testShape(self):
+        for pgp, gp in zip(self.pgps, self.gps):
+            self.assertAlmostEqual(
                 pgp.computeShape(),
                 gp.computeShape()
             )
 
+    def testResized(self):
+        for pgp, gp in zip(self.pgps, self.gps):
+            width, height = pgp.dimensions
             rpgp = pgp.resized(width+2, height+4)
+            # cppLib.resizedPsf calls Psf::resized, which redirects to
+            # PyGaussianPsf.resized above
+            rpgp2 = cppLib.resizedPsf(pgp, width+2, height+4)
             rgp = gp.resized(width+2, height+4)
             self.assertImagesAlmostEqual(
                 rpgp.computeImage(),
                 rgp.computeImage()
             )
-
-            cpgp = pgp.clone()
-            cgp = gp.clone()
-            self.assertIsNot(pgp, cpgp)
-            self.assertIsNot(gp, cgp)
-            self.assertImagesAlmostEqual(
-                cpgp.computeImage(),
-                cgp.computeImage()
-            )
             self.assertImagesAlmostEqual(
-                cpgp.computeImage(),
-                pgp.computeImage()
+                rpgp2.computeImage(),
+                rgp.computeImage()
             )
 
+    def testClone(self):
+        """Test different ways of invoking PyGaussianPsf.__deepcopy__
+        """
+        for pgp in self.pgps:
+            # directly
+            p1 = deepcopy(pgp)
+            # cppLib::clonedPsf -> Psf::clone
+            p2 = cppLib.clonedPsf(pgp)
+            # cppLib::clonedStorablePsf -> Psf::cloneStorable
+            p3 = cppLib.clonedStorablePsf(pgp)
+            # Psf::clone()
+            p4 = pgp.clone()
+
+            for p in [p1, p2, p3, p4]:
+                self.assertIsNot(pgp, p)
+                self.assertImagesEqual(
+                    pgp.computeImage(),
+                    p.computeImage()
+                )
+
+
+# Psf with position-dependent image, but nonetheless may use isFixed=True.
+# When isFixed=True, first image returned is cached for all subsequent image
+# queries
+class TestPsf(Psf):
+    def __init__(self, isFixed):
+        Psf.__init__(self, isFixed=isFixed)
+
+    def _doComputeKernelImage(self, position=None, color=None):
+        bbox = Box2I(Point2I(-3, -3), Extent2I(7, 7))
+        img = Image(bbox, dtype=np.float64)
+        x, y = np.ogrid[bbox.minY:bbox.maxY+1, bbox.minX:bbox.maxX+1]
+        rsqr = x**2 + y**2
+        if position.x >= 0.0:
+            img.array[:] = np.exp(-0.5*rsqr)
+        else:
+            img.array[:] = np.exp(-0.5*rsqr/4)
+        img.array /= np.sum(img.array)
+        return img
+
+
+class FixedPsfTestSuite(lsst.utils.tests.TestCase):
+    def setUp(self):
+        self.fixedPsf = TestPsf(isFixed=True)
+        self.floatPsf = TestPsf(isFixed=False)
+
+    def testFloat(self):
+        pos1 = Point2D(1.0, 1.0)
+        pos2 = Point2D(-1.0, -1.0)
+        img1 = self.floatPsf.computeKernelImage(pos1)
+        img2 = self.floatPsf.computeKernelImage(pos2)
+        self.assertFloatsNotEqual(img1.array, img2.array)
+
+    def testFixed(self):
+        pos1 = Point2D(1.0, 1.0)
+        pos2 = Point2D(-1.0, -1.0)
+        img1 = self.fixedPsf.computeKernelImage(pos1)
+        # Although _doComputeKernelImage would return a different image here due
+        # do the difference between pos1 and pos2, for the fixed Psf, the
+        # caching mechanism intercepts instead and _doComputeKernelImage is
+        # never called with position=pos2.  So img1 == img2.
+        img2 = self.fixedPsf.computeKernelImage(pos2)
+        self.assertFloatsEqual(img1.array, img2.array)
+
 
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass