Merge pull request #8 from lsst/tickets/DM-13065

DM-13065: Add low-level polynomials library.

.gitignore

@@ -21,6 +21,7 @@ doc/doxygen.conf
 examples/transform
 tests/.tests
 tests/.cache
+tests/test_polynomials
 tests/test_coordinates
 tests/test_spherePoint
 version.py

include/lsst/geom/polynomials.h

@@ -0,0 +1,103 @@
+// -*- LSST-C++ -*-
+/*
+ * 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 LSST_AFW_MATH_polynomials_h_INCLUDED
+#define LSST_AFW_MATH_polynomials_h_INCLUDED
+
+#include "lsst/geom/polynomials/BinomialMatrix.h"
+#include "lsst/geom/polynomials/PolynomialFunction1d.h"
+#include "lsst/geom/polynomials/Chebyshev1Function1d.h"
+#include "lsst/geom/polynomials/PolynomialFunction2d.h"
+#include "lsst/geom/polynomials/Chebyshev1Function2d.h"
+
+namespace lsst { namespace geom {
+
+/**
+ *  @namespace lsst::geom::polynomials Low-level polynomials (including special polynomials) in C++.
+ *
+ *  The geom::polynomials library provides low-level classes for
+ *  efficiently evaluating polynomial basis functions and expansions in C++.
+ *  The classes here:
+ *   - are not available in Python;
+ *   - are not polymorphic (no virtual functions);
+ *   - provide workspace objects to minimize memory allocations when
+ *     appropriate;
+ *   - do not throw exceptions (users are responsible for providing valid
+ *     inputs);
+ *   - have almost no outside dependencies (just Eigen);
+ *   - use templates to allow them to work with any array/vector objects (not
+ *     just Eigen).
+ *
+ *  They are intended to be used as the building blocks of higher-level
+ *  objects that are visible to Python users and persisted with our data
+ *  products, such as afw::math::ChebyshevBoundedField and the
+ *  afw::math::Function hierarchy.
+ *
+ *  At present, the library only includes support for 1-d and 2-d standard
+ *  polynomials and Chebyshev polynomials of the first kind, but adding
+ *  support for any other function defined by a recurrence relation (i.e. any
+ *  other special polynomial) should be extremely easy (see
+ *  RecurrenceBasis1d), and need not be done within the polynomials library
+ *  itself.
+ *
+ *  For both 1-d and 2-d, the library contains the following kinds of objects:
+ *
+ *   - Basis1d and Basis2d: objects that evaluate basis functions at
+ *     individual points.  The only concrete Basis1d implementation provided
+ *     at present is RecurrenceBasis1d template class, which can be used (with
+ *     the recurrence relation as a template parameter) to implement most
+ *     special polynomials.  Recurrence relations for standard polynomials
+ *     (PolynomialRecurrence) and Chebyshev polynomials of the first kind
+ *     (Chebyshev1Recurrence) are provided here.   The PackedBasis2d template
+ *     provides an implementation of Basis2d that evaluates a Basis1d over
+ *     each dimension.
+ *     @ref PolynomialBasis1d, @ref Chebyshev1Basis1d, @ref PolynomialBasis2d,
+ *     and @ref Chebyshev1Basis2d provide typedefs to instantiations of these
+ *     that should generally be preferred to explicit use of these templates.
+ *
+ *   - Function1d and Function2d: templates that combine a Basis1d or Basis2d
+ *     with a vector of coefficients.
+ *
+ *   - Scaling1d and Scaling2d: transformation objects that, with the
+ *     ScaledBasis1d and ScaledBasis2d templates, allow a basis or function to
+ *     be constructed that remaps points as part of evaluating the basis
+ *     functions.  For example, a @ref ScaledChebyshev1Basis1d combines both
+ *     the a Chebyshev polynomial basis and the scaling from some domain to
+ *     [-1, 1] that enables most special properties of Chebyshevs.  When
+ *     necessary, the simplified() functions can be used to convert a
+ *     ScaledPolynomialFunction1d or ScaledPolynomialFunction2d into an
+ *     equivalent PolynomialFunction1d or PolynomialFunction2d by folding the
+ *     scaling into the coefficients.
+ *
+ *  The library also includes a few utility classes and functions:
+ *
+ *   - BinomialMatrix provides an efficient and stable way to get
+ *     binomial coefficients.
+ *
+ *   - PackedIndexIterator and PackedIndexRange implement PackingOrders that
+ *     allow a conceptual 2-d triangular matrix to be implemented on top of
+ *     a 1-d array.
+ */
+namespace polynomials {
+
+}}} // namespace lsst::geom::polynomials
+
+#endif // !LSST_AFW_MATH_polynomials_h_INCLUDED

include/lsst/geom/polynomials/Basis1d.h

@@ -0,0 +1,108 @@
+// -*- LSST-C++ -*-
+/*
+ * 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 LSST_AFW_MATH_POLYNOMIALS_Basis1d_h_INCLUDED
+#define LSST_AFW_MATH_POLYNOMIALS_Basis1d_h_INCLUDED
+#ifdef DOXYGEN
+
+namespace lsst { namespace geom { namespace polynomials {
+
+/**
+ *  A basis interface for 1-d series expansions.
+ *
+ *  @note This class is only present in the documentation, as it represents an
+ *        abstract interface for which C++ (prior to C++20, at least) has no
+ *        language support.  It may be formalized into a true Concept when
+ *        that language feature is available.
+ */
+class Basis1d {
+public:
+
+    /// A Function1d object that uses this basis.
+    using Function = ...;
+
+    /// The type returned by scale().
+    using Scaled = ...;
+
+    /// Return the order of the basis.
+    std::size_t getOrder() const;
+
+    /// Return the number of elements in the basis.
+    std::size_t size() const;
+
+    /**
+     *  Return a scaled basis that delegates to a copy of `this`.
+     *
+     *  The scaled basis will transform all points by the given scaling
+     *  before evaluating the basis functions in the same way as `this`.
+     */
+    Scaled scaled(Scaling1d const & scaling) const;
+
+    /**
+     *  Evaluate a basis expansion with the given coefficients.
+     *
+     *  If the basis elements are @f$B_n(x)@f$ and the given coefficients are
+     *  a vector @f$a_n@f$, this computes
+     *  @f[
+     *      \sum_{n = 0}^{n \le N} a_n B_n(x)
+     *  @f]
+     *
+     *  @param[in] x             Point at which to evaluate the expansion.
+     *  @param[in] coefficients  Coefficients vector.  May be any type for
+     *                           which `coefficients[n]` returns an object
+     *                           convertible to `double` for all `n <=
+     *                           getOrder()`.  This includes
+     *                           `std::vector<double>`,
+     *                           `ndarray::Array<double,1>`,
+     *                           `Eigen::VectorXd`, and random access
+     *                           iterators.  If a lazy expression template
+     *                           object is passed, the elements of the
+     *                           expression will be evaluated only once.
+     *
+     *  @exceptsafe Does not throw unless `coefficients[n]` does, and provides
+     *              the same exception safety as it if it does.
+     */
+    template <typename Vector>
+    double sumWith(double x, Vector const & coefficients) const;
+
+    /**
+     *  Evaluate the basis at a given point.
+     *
+     *  @param[in] x       Point at which to evaluate the basis functions.
+     *  @param[out] basis  Output vector.  May be any type for which
+     *                     `coefficients[n]` returns a non-const reference to a
+     *                     floating-point value.  This includes
+     *                     `std::vector<double>`, `ndarray::Array<double,1>`,
+     *                     `Eigen::VectorXd`, `Eigen` view expressions, and
+     *                     mutable random access iterators.
+     *
+     *  @exceptsafe Does not throw unless `coefficients[n]` does, and provides
+     *              basic exception safety if it does.
+     */
+    template <typename Vector>
+    void fill(double x, Vector && basis) const;
+
+};
+
+}}} // namespace lsst::geom::polynomials
+
+#endif // DOXYGEN
+#endif // !LSST_AFW_MATH_POLYNOMIALS_Basis1d_h_INCLUDED

include/lsst/geom/polynomials/Basis2d.h

@@ -0,0 +1,107 @@
+// -*- LSST-C++ -*-
+/*
+ * 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 LSST_AFW_MATH_POLYNOMIALS_Basis2d_h_INCLUDED
+#define LSST_AFW_MATH_POLYNOMIALS_Basis2d_h_INCLUDED
+#ifdef DOXYGEN
+
+namespace lsst { namespace geom { namespace polynomials {
+
+/**
+ *  A basis interface for 2-d series expansions.
+ *
+ *  @note This class is only present in the documentation, as it represents an
+ *        abstract interface for which C++ (prior to C++20, at least) has no
+ *        language support.  It may be formalized into a true Concept when
+ *        that language feature is available.
+ */
+template <typename Basis1d>
+class Basis2d {
+public:
+
+    /// A Function2d object that uses this basis.
+    using Function = ...;
+
+    /// The type returned by scale().
+    using Scaled = ...;
+
+    /// The type returned by makeWorkspace().
+    using Workspace = ...;
+
+    /// Return the maximum order of the basis.
+    std::size_t getOrder() const;
+
+    /// Return the number of basis functions.
+    std::size_t size() const;
+
+    /**
+     *  Return a scaled basis that delegates to a copy of `this`.
+     *
+     *  The scaled basis will transform all points by the given scaling
+     *  before evaluating the basis functions in the same way as `this`.
+     */
+    Scaled scaled(Scaling2d const & first) const;
+
+    /// Allocate workspace that can be passed to sumWith() and fill() to avoid repeated memory allocations.
+    Workspace makeWorkspace() const;
+
+    /**
+     *  Evaluate a basis expansion with the given coefficients.
+     *
+     *  If the 1-d basis elements are @f$B_n(x)@f$ and the given coefficients are
+     *  a vector @f$a_{p, q}@f$, this computes
+     *  @f[
+     *      \sum_{p = 0, q = 0}^{p + q \le N} a_{p,q} B_{p}(x) B_{q}(y)
+     *  @f]
+     *
+     *  @param[in] point         Point at which to evaluate the expansion.
+     *  @param[in] coefficients  Flattened coefficients vector.
+     *                           See Basis1d::sumWith for more information.
+     */
+    template <typename Vector>
+    double sumWith(geom::Point2D const & point, Vector const & coefficients) const;
+
+    /// Evaluate a basis expansion with the given coefficients (external workspace version).
+    template <typename Vector>
+    double sumWith(geom::Point2D const & point, Vector const & coefficients, Workspace & workspace) const;
+
+    /**
+     *  Evaluate the basis at a given point.
+     *
+     *  @param[in] point   Point at which to evaluate the basis functions.
+     *  @param[out] basis  Flattened output vector.
+     *                     See Basis1d::fill more information.
+     */
+    template <typename Vector>
+    void fill(geom::Point2D const & point, Vector && basis) const;
+
+    /// Evaluate the basis at a given point (external workspace version).
+    template <typename Vector>
+    void fill(geom::Point2D const & point, Vector && basis, Workspace & workspace) const;
+
+private:
+    Basis1d _basis1d;
+};
+
+}}} // namespace lsst::geom::polynomials
+
+#endif // DOXYGEN
+#endif // !LSST_AFW_MATH_POLYNOMIALS_Basis2d_h_INCLUDED