RecurrenceBasis1d.h

Go to the documentation of this file.

// -*- 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_RecurrenceBasis1d_h_INCLUDED
#define LSST_AFW_MATH_POLYNOMIALS_RecurrenceBasis1d_h_INCLUDED

#include "lsst/geom/polynomials/SafeSum.h"
#include "lsst/geom/polynomials/ScaledBasis1d.h"

namespace lsst { namespace geom { namespace polynomials {

template <typename Basis>
class Function1d;

#ifdef DOXYGEN

struct Recurrence {

    static double getB0(double x);

    static double getB1(double x);

    static double next(double x, std::size_t n, double current, double previous);

};

#endif  // DOXYGEN


template <typename Recurrence>
class RecurrenceBasis1d {
public:

    using Function = Function1d<RecurrenceBasis1d>;

    using Scaled = ScaledBasis1d<RecurrenceBasis1d>;

    explicit RecurrenceBasis1d(std::size_t order) noexcept :
        _order(order)
    {}

    RecurrenceBasis1d(RecurrenceBasis1d const &) = default;

    RecurrenceBasis1d(RecurrenceBasis1d &&) = default;

    RecurrenceBasis1d & operator=(RecurrenceBasis1d const &) = default;

    RecurrenceBasis1d & operator=(RecurrenceBasis1d &&) = default;

    std::size_t getOrder() const noexcept { return _order; }

    std::size_t size() const noexcept { return _order + 1; }

    Scaled scaled(Scaling1d const & scaling) const noexcept {
        return Scaled(*this, scaling);
    }

    template <typename Vector>
    double sumWith(double x, Vector const & coefficients, SumMode mode=SumMode::FAST) const {
        // This universal lambda lets us effectively template most of the
        // implementation of this function on double vs. SafeSum<double>
        // without having to define an external template.
        auto accumulate = [x, coefficients, this](auto & sum) {
            double previous = Recurrence::getB0(x);
            if (_order > 0u) {
                double current = Recurrence::getB1(x);
                sum += coefficients[1]*current;
                for (std::size_t n = 2; n <= _order; ++n) {
                    double next = Recurrence::next(x, n, current, previous);
                    sum += coefficients[n]*next;
                    previous = current;
                    current = next;
                }
            }
        };
        double result = 0.0;
        if (mode == SumMode::FAST) {
            double z = Recurrence::getB0(x)*coefficients[0];
            accumulate(z);
            result = z;
        } else {
            SafeSum<double> z(Recurrence::getB0(x)*coefficients[0]);
            accumulate(z);
            result = static_cast<double>(z);
        }
        return result;
    }

    template <typename Vector>
    void fill(double x, Vector && basis) const {
        std::forward<Vector>(basis)[0] = Recurrence::getB0(x);
        if (_order > 0u) {
            std::forward<Vector>(basis)[1] = Recurrence::getB1(x);
            for (std::size_t n = 2; n <= _order; ++n) {
                std::forward<Vector>(basis)[n] = Recurrence::next(
                    x, n,
                    std::forward<Vector>(basis)[n - 1],
                    std::forward<Vector>(basis)[n - 2]
                );
            }
        }
    }

private:
    std::size_t _order;
};

}}} // namespace lsst::geom::polynomials

#endif // !LSST_AFW_MATH_POLYNOMIALS_RecurrenceBasis1d_h_INCLUDED