DM-52986: Add wrapper for SplineMap

include/astshim.h

@@ -41,6 +41,7 @@
 #include "astshim/ShiftMap.h"
 #include "astshim/SlaMap.h"
 #include "astshim/SphMap.h"
+#include "astshim/SplineMap.h"
 #include "astshim/TimeMap.h"
 #include "astshim/TranMap.h"
 #include "astshim/UnitMap.h"

include/astshim/SplineMap.h

@@ -0,0 +1,118 @@
+/*
+ * This file is part of astshim.
+ *
+ * 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/>.
+ */
+#ifndef ASTSHIM_SPLINEMAP_H
+#define ASTSHIM_SPLINEMAP_H
+
+#include <memory>
+#include <vector>
+
+#include "astshim/base.h"
+#include "astshim/Mapping.h"
+
+namespace ast {
+
+/**
+SplineMap is a @ref Mapping which performs a 2-in, 2-out Mapping defined by a pair of two-
+dimensional B-spline surfaces. The inverse transformation is fit iteratively, and requires that
+the output axes are a small (i.e. much less than the range of the input values) perturbation of the
+input axes. 
+### Attributes
+
+All those of @ref Mapping plus:
+
+- "InvNiter": Maximum number of iterations for iterative inverse.
+- "InvTol": Target relative error for iterative inverse.
+- "OutUnit": Determines how out-of-bounds input positions are handled. If 
+true, input points outside the bounding box are copied unchanged to output. If false, the output
+for these points will be nan.
+*/
+class SplineMap : public Mapping {
+    friend class Object;
+
+public:
+    /**
+    Construct a @ref SplineMap with the forward transform specified by the supplied b-splines and
+    the inverse transform fit iteratively.
+
+    @param[in] kx  Polynomial order in X input direction.
+    @param[in] ky  Polynomial order in Y input direction.
+    @param[in] nx  Number of coefficients in X input direction.
+    @param[in] ny  Number of coefficients in Y input direction
+    @param[in] tx  Array of knot positions for X input.
+    @param[in] ty  Array of knot positions for Y input.
+    @param[in] cu  Array of coefficients defining U output.
+    @param[in] cv  Array of coefficients defining V output.
+    @param[in] options Comma-separated list of attribute assignments.
+
+    @throws std::invalid_argument if the polynomial order of number of coefficients in either 
+        direction are negative.
+    @throws std::invalid_argument if the polynomial order plus the number of coefficients does
+        not equal the number of knot positions for a given direction.
+    @throws std::invalid_argument if the product of the number of coefficients in each direction
+        does not equal the number of coefficients provided.
+
+    The coefficients cu and cv are 1-dimensional arrays.
+    */
+    explicit SplineMap(int kx, int ky, int nx, int ny, std::vector<double> const &tx, std::vector<double> const &ty,
+                       std::vector<double> const &cu, std::vector<double> const &cv, std::string const &options = "")
+            : Mapping(reinterpret_cast<AstMapping *>(
+                      _makeRawSplineMap(kx, ky, nx, ny, tx, ty, cu, cv, options))) {}
+
+    ~SplineMap() override = default;
+
+    /// Copy constructor: make a deep copy
+    SplineMap(SplineMap const &) = default;
+    SplineMap(SplineMap &&) = default;
+    SplineMap &operator=(SplineMap const &) = delete;
+    SplineMap &operator=(SplineMap &&) = default;
+
+    /// Return a deep copy of this object.
+    std::shared_ptr<SplineMap> copy() const { return std::static_pointer_cast<SplineMap>(copyPolymorphic()); }
+
+    /// Get "InvNiter": Get maximum number of iterations for iterative inverse.
+    int getInvNiter() const { return getI("InvNiter"); }
+
+    /// Get "OutUnit": Is the out-of-bounds position transform an identity map?
+    int getOutUnit() const { return getI("OutUnit"); }
+
+    /// Get "InvTol": Get target relative error for iterative inverse.
+    double getInvTol() const { return getD("InvTol"); }
+
+
+protected:
+    std::shared_ptr<Object> copyPolymorphic() const override {
+        return copyImpl<SplineMap, AstSplineMap>();
+    }
+
+    /// Construct a SplineMap from an raw AST pointer
+    explicit SplineMap(AstSplineMap *map);
+
+private:
+    /// Make a raw AstSplineMap with a specified forward transform and an iterative inverse.
+    AstSplineMap *_makeRawSplineMap(int kx, int ky, int nx, int ny, std::vector<double> const &tx, std::vector<double> const &ty,
+                       std::vector<double> const &cu, std::vector<double> const &cv, std::string const &options = "") const;
+};
+
+}  // namespace ast
+
+#endif

python/astshim/SConscript

@@ -33,6 +33,7 @@ scripts.BasicSConscript.python(['_astshimLib'], [
     "slaMap.cc",
     "specFrame.cc",
     "sphMap.cc",
+    "splineMap.cc",
     "stream.cc",
     "table.cc",
     "timeFrame.cc",

python/astshim/_astshimLib.cc

@@ -61,6 +61,7 @@ void wrapSkyFrame(WrapperCollection&);
 void wrapSlaMap(WrapperCollection&);
 void wrapSpecFrame(WrapperCollection&);
 void wrapSphMap(WrapperCollection&);
+void wrapSplineMap(WrapperCollection&);
 void wrapStream(WrapperCollection&);
 void wrapTable(WrapperCollection&);
 void wrapTimeFrame(WrapperCollection&);
@@ -110,6 +111,7 @@ PYBIND11_MODULE(_astshimLib, mod) {
     wrapSlaMap(wrappers);
     wrapSpecFrame(wrappers);
     wrapSphMap(wrappers);
+    wrapSplineMap(wrappers);
     wrapStream(wrappers);
     wrapTimeFrame(wrappers);
     wrapTimeMap(wrappers);

python/astshim/splineMap.cc

@@ -0,0 +1,52 @@
+/*
+ * This file is part of astshim.
+ *
+ * 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/>.
+ */
+#include <memory>
+#include <vector>
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "lsst/cpputils/python.h"
+
+#include "astshim/Mapping.h"
+#include "astshim/SplineMap.h"
+
+namespace py = pybind11;
+using namespace pybind11::literals;
+
+namespace ast {
+void wrapSplineMap(lsst::cpputils::python::WrapperCollection &wrappers) {
+    using PySplineMap=py::classh<SplineMap, Mapping>;
+    wrappers.wrapType(PySplineMap(wrappers.module, "SplineMap"), [](auto &mod, auto &cls) {
+
+        cls.def(py::init<int, int, int, int, std::vector<double> const &, std::vector<double> const &,
+                       std::vector<double> const &, std::vector<double> const &, std::string const &>(), "kx"_a, 
+                       "ky"_a, "nx"_a, "ny"_a, "tx"_a, "ty"_a, "cu"_a, "cv"_a, "options"_a = "");
+        cls.def(py::init<SplineMap const &>());
+        cls.def("copy", &SplineMap::copy);
+        cls.def_property_readonly("invNIter", &SplineMap::getInvNiter);
+        cls.def_property_readonly("outUnit", &SplineMap::getOutUnit);
+        cls.def_property_readonly("invTol", &SplineMap::getInvTol);
+    });
+}
+
+}  // namespace ast

src/Object.cc

@@ -54,6 +54,7 @@
 #include "astshim/SlaMap.h"
 #include "astshim/SpecFrame.h"
 #include "astshim/SphMap.h"
+#include "astshim/SplineMap.h"
 #include "astshim/Table.h"
 #include "astshim/TimeFrame.h"
 #include "astshim/TimeMap.h"
@@ -114,6 +115,7 @@ std::shared_ptr<Object> Object::_basicFromAstObject(AstObject *rawObj) {
                     {"SlaMap", makeShim<SlaMap, AstSlaMap>},
                     {"SpecFrame", makeShim<SpecFrame, AstSpecFrame>},
                     {"SphMap", makeShim<SphMap, AstSphMap>},
+                    {"SplineMap", makeShim<SplineMap, AstSplineMap>},
                     {"TimeFrame", makeShim<TimeFrame, AstTimeFrame>},
                     {"Table", makeShim<Table, AstTable>},
                     {"TimeMap", makeShim<TimeMap, AstTimeMap>},

src/SplineMap.cc

@@ -0,0 +1,64 @@
+/*
+ * This file is part of astshim.
+ *
+ * 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/>.
+ */
+#include <sstream>
+#include <stdexcept>
+
+#include "astshim/SplineMap.h"
+
+namespace ast {
+
+SplineMap::SplineMap(AstSplineMap *map) : Mapping(reinterpret_cast<AstMapping *>(map)) {
+    if (!astIsASplineMap(getRawPtr())) {
+        std::ostringstream os;
+        os << "this is a " << getClassName() << ", which is not a SplineMap";
+        throw std::invalid_argument(os.str());
+    }
+}
+
+/// Make a raw AstSplineMap with a specified forward transform.
+AstSplineMap *SplineMap::_makeRawSplineMap(int kx, int ky, int nx, int ny, std::vector<double> const &tx,
+    std::vector<double> const &ty, std::vector<double> const &cu, std::vector<double> const &cv,
+    std::string const &options) const {
+
+    const size_t nTx = tx.size();
+    const size_t nTy = ty.size();
+    const size_t nCoeffsU = cu.size();
+    const size_t nCoeffsV = cv.size();
+
+    if ((kx < 0) || (ky < 0) || (nx < 0) || (ny < 0)) {
+        throw std::invalid_argument("The polynomial order and number of coefficients must not be negative.");
+    }
+    if (((size_t)(kx + nx) != nTx) || ((size_t)(ky + ny) != nTy)) {
+        throw std::invalid_argument("The length of the knot positions must equal the polynomial order plus "
+            "the number of coefficients.");
+    }
+    if (((size_t)(nx * ny) != nCoeffsU) || ((size_t)(nx * ny) != nCoeffsV)) {
+        throw std::invalid_argument("The length of the coefficients must equal product of the arguments nx "
+            "and ny.");
+    }
+
+    return reinterpret_cast<AstSplineMap *>(astSplineMap(kx, ky, nx, ny, tx.data(), ty.data(), cu.data(),
+        cv.data(), options.c_str()));
+}
+
+}  // namespace ast

tests/test_splineMap.py

@@ -0,0 +1,99 @@
+import unittest
+
+import numpy as np
+import numpy.testing as npt
+from scipy.interpolate import RectBivariateSpline
+
+import astshim as ast
+from astshim.test import MappingTestCase
+
+
+class TestSplineMap(MappingTestCase):
+
+    def setUp(self):
+        self.nx = 150
+        self.ny = 100
+        self.k = 4
+        a = 4.0
+        b = 0.2
+        c = 3.0
+
+        self.xs = np.linspace(1, self.nx, self.nx)
+        self.ys = np.linspace(1, self.ny, self.ny)
+
+        fcx = np.zeros((self.nx, self.ny))
+        fcy = np.zeros((self.nx, self.ny))
+        for j in range(self.ny):
+            y = self.ys[j]
+            for i in range(self.nx):
+                x = self.xs[i]
+                fcx[i][j] = x + a * np.sin(b * x) * np.cos(b * y)
+                fcy[i][j] = y + c * np.cos(b * x) * np.sin(b * y)
+
+        self.splinex = RectBivariateSpline(self.xs, self.ys, fcx, s=0)
+        (self.tx, self.ty) = self.splinex.get_knots()
+        self.arx = self.splinex.get_coeffs()
+
+        self.spliney = RectBivariateSpline(self.xs, self.ys, fcy, s=0)
+        self.ary = self.spliney.get_coeffs()
+
+        self.splineMap = ast.SplineMap(
+            self.k, self.k, self.nx, self.ny, self.tx, self.ty, self.arx, self.ary
+        )
+
+    def test_SplineMap(self):
+        """Test that the forward and inverse transforms match
+        scipy.interpolate.RectBivariateSpline.
+        """
+        xval = np.array([0.0, 12.5, 12.0, 1.0, 15.0, 1.0, 95.0, 149.5, 151.2, 77.77])
+        yval = np.array([-1.0, 8.8, 8.0, 1.0, 15.0, 76.0, 100.0, 99.8, 82.3, 54.3])
+
+        u = self.splinex.ev(xval, yval)
+        v = self.spliney.ev(xval, yval)
+
+        outData = self.splineMap.applyForward(np.array([xval, yval]))
+
+        outOfBounds = (
+            (xval > self.xs.max())
+            | (xval < self.xs.min())
+            | (yval < self.ys.min())
+            | (yval > self.ys.max())
+        )
+
+        npt.assert_equal(outData[:, outOfBounds], np.nan)
+        npt.assert_almost_equal(outData[0, ~outOfBounds], u[~outOfBounds])
+        npt.assert_almost_equal(outData[1, ~outOfBounds], v[~outOfBounds])
+
+        reverseTrip = self.splineMap.applyInverse(outData[:, ~outOfBounds].copy())
+
+        npt.assert_almost_equal(reverseTrip[0], xval[~outOfBounds])
+        npt.assert_almost_equal(reverseTrip[1], yval[~outOfBounds])
+
+    def test_splineMapAttributes(self):
+        """Check that SplineMap attributes can be accessed and return expected
+        values."""
+        self.assertEqual(self.splineMap.invTol, 1e-6)
+        self.assertEqual(self.splineMap.invNIter, 6)
+        self.assertEqual(self.splineMap.outUnit, False)
+
+        # Initialize SplineMap with other attributes and check they are set
+        # correctly.
+        options = "InvNIter=5,InvTol=1e-7,OutUnit=1"
+        splineMap = ast.SplineMap(
+            self.k,
+            self.k,
+            self.nx,
+            self.ny,
+            self.tx,
+            self.ty,
+            self.arx,
+            self.ary,
+            options=options,
+        )
+        self.assertEqual(splineMap.invTol, 1e-7)
+        self.assertEqual(splineMap.invNIter, 5)
+        self.assertEqual(splineMap.outUnit, 1)
+
+
+if __name__ == "__main__":
+    unittest.main()