DM-42920: Add vectorized evaluateArray to CoaddBoundedField.

include/lsst/meas/algorithms/CoaddBoundedField.h

@@ -80,6 +80,9 @@ class CoaddBoundedField : public afw::table::io::PersistableFacade<CoaddBoundedF
     /// Get the elements vector
     ElementVector getElements() const { return _elements; };
 
+    /// Get the throwOnMissing setting
+    bool getThrowOnMissing() const { return _throwOnMissing; };
+
     /**
      *  @brief Return true if the CoaddBoundedField persistable (always true).
      */

python/lsst/meas/algorithms/__init__.py

@@ -63,6 +63,7 @@
 from .noise_covariance import *
 from .reinterpolate_pixels import *
 from .setPrimaryFlags import *
+from .coaddBoundedField import *
 
 from .version import *
 

python/lsst/meas/algorithms/coaddBoundedField.cc

@@ -75,10 +75,11 @@ void declareCoaddBoundedField(lsst::cpputils::python::WrapperCollection &wrapper
         cls.def("__imul__", &CoaddBoundedField::operator*);
 
         /* Members */
-        cls.def("evaluate", &CoaddBoundedField::evaluate);
+        cls.def("_evaluate", &CoaddBoundedField::evaluate);
         cls.def("getCoaddWcs", &CoaddBoundedField::getCoaddWcs);
         cls.def("getDefault", &CoaddBoundedField::getDefault);
         cls.def("getElements", &CoaddBoundedField::getElements);
+        cls.def("getThrowOnMissing", &CoaddBoundedField::getThrowOnMissing);
     });
     afw::table::io::python::addPersistableMethods<CoaddBoundedField>(clsField);
 }

python/lsst/meas/algorithms/coaddBoundedField.py

@@ -0,0 +1,86 @@
+# This file is part of meas_algorithms.
+#
+# 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/>.
+
+import numpy as np
+
+from lsst.geom import Point2D
+from lsst.utils import continueClass
+
+from ._algorithmsLib import CoaddBoundedField
+
+__all__ = ["CoaddBoundedField"]
+
+
+@continueClass
+class CoaddBoundedField:  # noqa: F811
+    def evaluate(self, x, y=None):
+        """Evaluate the CoaddBoundedField.
+
+        This accepts either a Point2D or an array of x and y
+        positions.
+
+        When arrays are passed, this uses a vectorized version
+        of CoaddBoundedField::evaluate(). If the coadd bounded
+        field has throwOnMissing then this will return NaN
+        for missing values; otherwise it will return 0.0.
+
+        Parameters
+        ----------
+        x : `lsst.geom.Point2D` or `np.ndarray`
+            Array of x values.
+        y : `np.ndarray`, optional
+            Array of y values.
+
+        Returns
+        -------
+        values : `float` or `np.ndarray`
+            Evaluated value or array of values.
+        """
+        if isinstance(x, Point2D):
+            return self._evaluate(x)
+
+        _x = np.atleast_1d(x).ravel()
+        _y = np.atleast_1d(y).ravel()
+
+        if len(_x) != len(_y):
+            raise ValueError("x and y arrays must be the same length.")
+
+        ra, dec = self.getCoaddWcs().pixelToSkyArray(_x, _y)
+
+        sums = np.zeros(len(_x))
+        wts = np.zeros_like(sums)
+
+        for elt in self.getElements():
+            ix, iy = elt.wcs.skyToPixelArray(ra, dec)
+            in_box = elt.field.getBBox().contains(ix, iy)
+            if elt.validPolygon:
+                in_box &= elt.validPolygon.contains(ix, iy)
+            sums[in_box] += elt.weight*elt.field.evaluate(ix[in_box], iy[in_box])
+            wts[in_box] += elt.weight
+
+        good = (wts > 0)
+        values = np.zeros(len(_x))
+        values[good] = sums[good]/wts[good]
+
+        if self.getThrowOnMissing():
+            values[~good] = np.nan
+
+        return values

tests/test_coaddApCorrMap.py

@@ -114,20 +114,48 @@ def testCoaddApCorrMap(self):
         os.unlink(filename)
 
     def assertApCorrMap(self, apCorrMap, pointList):
+        # Check point-by-point.
+        pointsX = np.zeros(len(pointList))
+        pointsY = np.zeros(len(pointList))
+        expectedArray = np.zeros(len(pointList))
         for i, point in enumerate(pointList):
             weights = [i + 1, i + 2]
             values = [i + 1, i + 2]
             expected = sum((w*v for w, v in zip(weights, values)), 0.0)/sum(weights)
             actual = apCorrMap["only"].evaluate(point)
             self.assertEqual(actual, expected)
+            # Create arrays for vectorized call.
+            expectedArray[i] = expected
+            pointsX[i] = point.getX()
+            pointsY[i] = point.getY()
+
+        # Check vectorized
+        actualArray = apCorrMap["only"].evaluate(pointsX, pointsY)
+        self.assertFloatsAlmostEqual(actualArray, expectedArray)
+
+        # Check vectorized single point.
+        actualArray0 = apCorrMap["only"].evaluate(pointsX[0], pointsY[0])
+        self.assertFloatsAlmostEqual(actualArray0, expectedArray[0])
 
     def assertApCorrMapValid(self, apCorrMap, pointList):
+        # Check point-by-point.
+        pointsX = np.zeros(len(pointList))
+        pointsY = np.zeros(len(pointList))
+        expectedArray = np.zeros(len(pointList))
         for i, point in enumerate(pointList):
             weights = [i + 2]
             values = [i + 2]
             expected = sum((w*v for w, v in zip(weights, values)), 0.0)/sum(weights)
             actual = apCorrMap["only"].evaluate(point)
             self.assertEqual(actual, expected)
+            # Create arrays for vectorized call.
+            expectedArray[i] = expected
+            pointsX[i] = point.getX()
+            pointsY[i] = point.getY()
+
+        # Check vectorized
+        actualArray = apCorrMap["only"].evaluate(pointsX, pointsY)
+        self.assertFloatsAlmostEqual(actualArray, expectedArray)
 
 
 class TestMemory(lsst.utils.tests.MemoryTestCase):