Add Transform::getJacobian.

.gitignore

@@ -136,6 +136,7 @@ tests/statistics
 tests/statisticsSpeed
 tests/testLib.py
 tests/testSpherePoint
+tests/test_transform
 tests/testTrapezoidalPacker
 tests/testWcs
 tests/testWarpGpu

include/lsst/afw/geom/Transform.h

@@ -28,6 +28,7 @@
 #include <vector>
 
 #include "astshim.h"
+#include "Eigen/Core"
 #include "ndarray.h"
 
 #include "lsst/afw/geom/Endpoint.h"
@@ -167,6 +168,27 @@ class Transform {
      */
     Transform<ToEndpoint, FromEndpoint> getInverse() const;
 
+    /**
+     * The Jacobian matrix of this Transform.
+     *
+     * The matrix is defined only if this object has a forward transform.
+     *
+     * @param x the position at which the Jacobian shall be evaluated
+     * @returns a matrix `J` with `getToEndpoint().getNAxes()` rows and
+     *          `getFromEndpoint().getNAxes()` columns. `J(i,j)` shall be the
+     *          rate of change of the `i`th output coordinate with respect to
+     *          the `j`th input coordinate, or `NaN` if the derivative cannot
+     *          be calculated.
+     *
+     * @exceptsafe Provides basic exception safety.
+     *
+     * @note The derivatives may be estimated by sampling and interpolating
+     *       this Transform in the neighborhood of `x`. If the implementation
+     *       requires interpolation, computation of the Jacobian may require
+     *       hundreds of evaluations of @ref tranForward.
+     */
+    Eigen::MatrixXd getJacobian(FromPoint const &x) const;
+
 private:
     FromEndpoint const _fromEndpoint;
     std::shared_ptr<const ast::FrameSet> _frameSet;

python/lsst/afw/geom/SConscript

@@ -13,7 +13,7 @@ scripts.BasicSConscript.pybind11([
     'linearTransform',
     'xyTransform',
 	'endpoint',
-	'transform',
+	'transform/transform',
     ],
     addUnderscore=False
 )

python/lsst/afw/geom/transform/__init__.py

@@ -0,0 +1,26 @@
+#
+# LSST Data Management System
+# Copyright 2008-2017 LSST/AURA.
+#
+# This product includes software developed by the
+# LSST Project (http://www.lsst.org/).
+#
+# 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 LSST License Statement and
+# the GNU General Public License along with this program.  If not,
+# see <http://www.lsstcorp.org/LegalNotices/>.
+#
+
+from __future__ import absolute_import
+
+from .transform import *
+from .transformContinued import *

python/lsst/afw/geom/transform.cc → python/lsst/afw/geom/transform/transform.cc

@@ -81,6 +81,7 @@ void declareTransform(py::module &mod, std::string const &fromName, std::string
     cls.def("tranInverse", (FromArray (Class::*)(ToArray const &) const) & Class::tranInverse, "array"_a);
     cls.def("tranInverse", (FromPoint (Class::*)(ToPoint const &) const) & Class::tranInverse, "point"_a);
     cls.def("getInverse", &Class::getInverse);
+    cls.def("_getJacobian", &Class::getJacobian);
     // str(self) = "<Python class name>[<nIn>-><nOut>]"
     cls.def("__str__", [pyClassName](Class const &self) { return formatStr(self, pyClassName); });
     // repr(self) = "lsst.afw.geom.<Python class name>[<nIn>-><nOut>]"

python/lsst/afw/geom/transform/transformContinued.py

@@ -0,0 +1,42 @@
+#
+# LSST Data Management System
+# Copyright 2017 LSST/AURA.
+#
+# This product includes software developed by the
+# LSST Project (http://www.lsst.org/).
+#
+# 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 LSST License Statement and
+# the GNU General Public License along with this program.  If not,
+# see <http://www.lsstcorp.org/LegalNotices/>.
+#
+from __future__ import absolute_import, division, print_function
+
+import sys
+
+# Too many classes for individually updating them to be practical
+from .transform import *
+
+__all__ = []
+
+
+def getJacobian(self, x):
+    # Force 2D matrix over numpy's protests
+    matrix = self._getJacobian(x)
+    matrix.shape = (self.getToEndpoint().getNAxes(),
+                    self.getFromEndpoint().getNAxes())
+    return matrix
+
+for name in dir():
+    if name.startswith("Transform"):
+        cls = getattr(sys.modules[__name__], name)
+        cls.getJacobian = getJacobian

src/geom/Transform.cc

@@ -20,6 +20,7 @@
  * see <http://www.lsstcorp.org/LegalNotices/>.
  */
 
+#include <exception>
 #include <memory>
 #include <ostream>
 #include <sstream>
@@ -30,6 +31,7 @@
 #include "lsst/afw/geom/Point.h"
 #include "lsst/afw/geom/SpherePoint.h"
 #include "lsst/afw/geom/Transform.h"
+#include "lsst/pex/exceptions/Exception.h"
 
 namespace lsst {
 namespace afw {
@@ -112,11 +114,30 @@ Transform<ToEndpoint, FromEndpoint> Transform<FromEndpoint, ToEndpoint>::getInve
         // don't throw std::bad_cast because it doesn't let you provide debugging info
         std::ostringstream buffer;
         buffer << "FrameSet.getInverse() does not return a FrameSet. Called from: " << _frameSet;
-        throw std::logic_error(buffer.str());
+        throw LSST_EXCEPT(pex::exceptions::LogicError, buffer.str());
     }
     return Transform<ToEndpoint, FromEndpoint>(*inverse);
 }
 
+template <typename FromEndpoint, typename ToEndpoint>
+Eigen::MatrixXd Transform<FromEndpoint, ToEndpoint>::getJacobian(FromPoint const &x) const {
+    try {
+        int const nIn = _fromEndpoint.getNAxes();
+        int const nOut = _toEndpoint.getNAxes();
+        std::vector<double> const point = _fromEndpoint.dataFromPoint(x);
+
+        Eigen::MatrixXd jacobian(nOut, nIn);
+        for (int i = 0; i < nOut; ++i) {
+            for (int j = 0; j < nIn; ++j) {
+                jacobian(i, j) = _frameSet->rate(point, i + 1, j + 1);
+            }
+        }
+        return jacobian;
+    } catch (std::bad_alloc const &e) {
+        std::throw_with_nested(LSST_EXCEPT(pex::exceptions::MemoryError, "Could not allocate Jacobian."));
+    }
+}
+
 template <typename FromEndpoint, typename ToEndpoint>
 std::ostream &operator<<(std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform) {
     auto const frameSet = transform.getFrameSet();

tests/test_transform.cc

@@ -0,0 +1,110 @@
+// -*- LSST-C++ -*-
+
+/*
+ * LSST Data Management System
+ * See COPYRIGHT file at the top of the source tree.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * 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 LSST License Statement and
+ * the GNU General Public License along with this program. If not,
+ * see <http://www.lsstcorp.org/LegalNotices/>.
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE TransformCpp
+
+#include <array>
+
+#include "boost/test/unit_test.hpp"
+
+#include "lsst/afw/geom/Transform.h"
+
+/*
+ * Unit tests for C++-only functionality in Transform.
+ *
+ * See test_transform.py for remaining unit tests.
+ */
+namespace lsst {
+namespace afw {
+namespace geom {
+
+/**
+ * Make an ast::PolyMap suitable for testing.
+ *
+ * The forward transform is as follows:
+ * @f[ f_j(\vec{x}) = \sum_{i} C_{ij} x_{i}^2 @f]
+ * where @f$ C_{ij} = 0.001 (i+j+1) @f$.
+ *
+ * The equation is chosen for the following reasons:
+ * - It is well defined for any value of `nIn`, `nOut`.
+ * - It stays small for small `x`, to avoid wraparound of angles for
+ *      SpherePoint endpoints.
+ *
+ * @param nIn,nOut The input and output dimensions of the desired PolyMap.
+ * @returns a Mapping with a forward transform described by the equation above
+ *          and no inverse transform.
+ */
+ast::PolyMap makeForwardPolyMap(size_t nIn, size_t nOut) {
+    using namespace ndarray;
+
+    double const baseCoeff = 0.001;
+    Array<double, 2, 2> forwardCoeffs = allocate(ndarray::makeVector(nOut * nIn, 2 + nIn));
+    for (size_t iOut = 0; iOut < nOut; ++iOut) {
+        double const coeffOffset = baseCoeff * iOut;
+        for (size_t iIn = 0; iIn < nIn; ++iIn) {
+            auto row = forwardCoeffs[iOut * nIn + iIn];
+            row[0] = baseCoeff * (iIn + 1) + coeffOffset;  // Coefficient
+            row[1] = iOut + 1;                             // Compute f_iOut
+            row[view(2, 2 + nIn)] = 0;                     // Ignore most variables
+            row[2 + iIn] = 2;                              // Square x_iIn
+        }
+    }
+
+    return ast::PolyMap(forwardCoeffs, nOut, "IterInverse=0");
+}
+
+/**
+ * Tests whether the result of SpherePoint::getJacobian(FromPoint const&)
+ * has the specified dimensions.
+ *
+ * The Python version of this method follows a slightly different spec to
+ * conform to the numpy convention that length=1 dimensions do not exist.
+ */
+BOOST_AUTO_TEST_CASE(getJacobianDimensions) {
+    using GenericTransform = Transform<GenericEndpoint, GenericEndpoint>;
+    std::array<size_t, 6> const dimensions = {1, 2, 3, 4, 5, 6};
+
+    for (auto nIn : dimensions) {
+        for (auto nOut : dimensions) {
+            std::string msg = " [nIn=" + std::to_string(nIn) + ", nOut=" + std::to_string(nOut) + "]";
+            auto polyMap = makeForwardPolyMap(nIn, nOut);
+            GenericTransform transform(polyMap);
+
+            // Don't care about elements, so zero initialization is ok
+            auto inPoint = std::vector<double>(nIn);
+            Eigen::MatrixXd jacobian = transform.getJacobian(inPoint);
+
+            auto fromAxes = transform.getFromEndpoint().getNAxes();
+            auto toAxes = transform.getToEndpoint().getNAxes();
+            BOOST_TEST(jacobian.rows() == toAxes, "Matrix has " << jacobian.rows() << " rows, expected "
+                                                                << toAxes << msg);
+            BOOST_TEST(jacobian.cols() == fromAxes, "Matrix has " << jacobian.cols() << " columns, expected "
+                                                                  << fromAxes << msg);
+        }
+    }
+}
+}
+}
+} /* namespace lsst::afw::geom */