PackedBasis2d.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_PackedBasis2d_h_INCLUDED
#define LSST_AFW_MATH_POLYNOMIALS_PackedBasis2d_h_INCLUDED

#include "lsst/geom/Point.h"
#include "lsst/geom/polynomials/SafeSum.h"
#include "lsst/geom/polynomials/PackedIndex.h"
#include "lsst/geom/polynomials/ScaledBasis2d.h"

namespace lsst { namespace geom { namespace polynomials {

template <typename Basis1d, PackingOrder packing>
class PackedBasis2d;


class PackedBasisWorkspace2d {
public:

    explicit PackedBasisWorkspace2d(std::size_t order) : _x(order + 1), _y(order + 1) {}

    std::size_t getOrder() const { return _x.size() - 1; }

private:

    template <typename Recurrence, PackingOrder packing>
    friend class PackedBasis2d;

    Eigen::VectorXd _x;
    Eigen::VectorXd _y;
};

template <typename Basis>
class Function2d;

template <typename Basis1d, PackingOrder packing>
class PackedBasis2d {
public:

    using Function = Function2d<PackedBasis2d>;

    using Scaled = ScaledBasis2d<PackedBasis2d>;

    using Workspace = PackedBasisWorkspace2d;

    using IndexRange = PackedIndexRange<packing>;

    static constexpr std::size_t computeSize(std::size_t order) { return IndexRange::computeSize(order); }

    explicit PackedBasis2d(Basis1d const & basis1d) : _basis1d(basis1d) {}

    template <typename ...Args>
    explicit PackedBasis2d(Args&& ...args) : _basis1d(std::forward<Args>(args)...) {}

    PackedBasis2d(PackedBasis2d const &) = default;

    PackedBasis2d(PackedBasis2d &&) = default;

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

    PackedBasis2d & operator=(PackedBasis2d &&) = default;

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

    std::size_t size() const noexcept{ return IndexRange::computeSize(getOrder()); }

    Scaled scaled(Scaling2d const & first) const {
        return Scaled(*this, first);
    }

    std::size_t index(std::size_t x, std::size_t y) const {
        return IndexRange::computeIndex(x, y);
    }

    IndexRange getIndices() const noexcept {
        return IndexRange(typename IndexRange::iterator(), IndexRange::iterator::makeEnd(getOrder()));
    }

    Workspace makeWorkspace() const { return Workspace(getOrder());}

    template <typename Vector>
    double sumWith(geom::Point2D const & point, Vector const & coefficients,
                   Workspace & workspace, SumMode mode=SumMode::FAST) const {
        assert(workspace.getOrder() >= getOrder());
        _basis1d.fill(point.getX(), workspace._x);
        _basis1d.fill(point.getY(), workspace._y);
        // 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 = [coefficients, &workspace, this](auto & sum) {
            for (auto const & index : getIndices()) {
                sum += coefficients[index.flat]*workspace._x[index.nx]*workspace._y[index.ny];
            }
        };
        double result = 0.0;
        if (mode == SumMode::FAST) {
            double z = 0.0;
            accumulate(z);
            result = z;
        } else {
            SafeSum<double> z;
            accumulate(z);
            result = static_cast<double>(z);
        }
        return result;
    }

    template <typename Vector>
    double sumWith(geom::Point2D const & point, Vector const & coefficients,
                   SumMode mode=SumMode::FAST) const {
        auto workspace = makeWorkspace();
        return sumWith(point, coefficients, workspace, mode);
    }

    template <typename Vector>
    void fill(geom::Point2D const & point, Vector && basis, Workspace & workspace) const {
        assert(workspace.getOrder() >= getOrder());
        _basis1d.fill(point.getX(), workspace._x);
        _basis1d.fill(point.getY(), workspace._y);
        for (auto const & index : getIndices()) {
            std::forward<Vector>(basis)[index.flat] = workspace._x[index.nx]*workspace._y[index.ny];
        }
    }

    template <typename Vector>
    void fill(geom::Point2D const & point, Vector && basis) const {
        auto workspace = makeWorkspace();
        fill(point, std::forward<Vector>(basis), workspace);
    }

private:
    Basis1d _basis1d;
};

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

#endif // !LSST_AFW_MATH_POLYNOMIALS_PackedBasis2d_h_INCLUDED