Loading [MathJax]/extensions/tex2jax.js
LSST Applications g0fba68d861+83433b07ee,g16d25e1f1b+23bc9e47ac,g1ec0fe41b4+3ea9d11450,g1fd858c14a+9be2b0f3b9,g2440f9efcc+8c5ae1fdc5,g35bb328faa+8c5ae1fdc5,g4a4af6cd76+d25431c27e,g4d2262a081+c74e83464e,g53246c7159+8c5ae1fdc5,g55585698de+1e04e59700,g56a49b3a55+92a7603e7a,g60b5630c4e+1e04e59700,g67b6fd64d1+3fc8cb0b9e,g78460c75b0+7e33a9eb6d,g786e29fd12+668abc6043,g8352419a5c+8c5ae1fdc5,g8852436030+60e38ee5ff,g89139ef638+3fc8cb0b9e,g94187f82dc+1e04e59700,g989de1cb63+3fc8cb0b9e,g9d31334357+1e04e59700,g9f33ca652e+0a83e03614,gabe3b4be73+8856018cbb,gabf8522325+977d9fabaf,gb1101e3267+8b4b9c8ed7,gb89ab40317+3fc8cb0b9e,gc0af124501+57ccba3ad1,gcf25f946ba+60e38ee5ff,gd6cbbdb0b4+1cc2750d2e,gd794735e4e+7be992507c,gdb1c4ca869+be65c9c1d7,gde0f65d7ad+c7f52e58fe,ge278dab8ac+6b863515ed,ge410e46f29+3fc8cb0b9e,gf35d7ec915+97dd712d81,gf5e32f922b+8c5ae1fdc5,gf618743f1b+747388abfa,gf67bdafdda+3fc8cb0b9e,w.2025.18
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
PackedBasis2d.h
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * Developed for the LSST Data Management System.
4 * This product includes software developed by the LSST Project
5 * (https://www.lsst.org).
6 * See the COPYRIGHT file at the top-level directory of this distribution
7 * for details of code ownership.
8 *
9 * This program is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 3 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program. If not, see <https://www.gnu.org/licenses/>.
21 */
22#ifndef LSST_AFW_MATH_POLYNOMIALS_PackedBasis2d_h_INCLUDED
23#define LSST_AFW_MATH_POLYNOMIALS_PackedBasis2d_h_INCLUDED
24
25#include "lsst/geom/Point.h"
29
30namespace lsst { namespace geom { namespace polynomials {
31
32template <typename Basis1d, PackingOrder packing>
33class PackedBasis2d;
34
35
41public:
42
44 explicit PackedBasisWorkspace2d(std::size_t order) : _x(order + 1), _y(order + 1) {}
45
47 std::size_t getOrder() const { return _x.size() - 1; }
48
49private:
50
51 template <typename Recurrence, PackingOrder packing>
52 friend class PackedBasis2d;
53
54 Eigen::VectorXd _x;
55 Eigen::VectorXd _y;
56};
57
58template <typename Basis>
59class Function2d;
60
74template <typename Basis1d, PackingOrder packing>
76public:
77
80
83
86
89
91 static constexpr std::size_t computeSize(std::size_t order) { return IndexRange::computeSize(order); }
92
94 explicit PackedBasis2d(Basis1d const & basis1d) : _basis1d(basis1d) {}
95
97 template <typename ...Args>
98 explicit PackedBasis2d(Args&& ...args) : _basis1d(std::forward<Args>(args)...) {}
99
101 PackedBasis2d(PackedBasis2d const &) = default;
102
105
108
111
113 std::size_t getOrder() const noexcept { return _basis1d.getOrder(); }
114
116 std::size_t size() const noexcept{ return IndexRange::computeSize(getOrder()); }
117
124 Scaled scaled(Scaling2d const & first) const {
125 return Scaled(*this, first);
126 }
127
130 return IndexRange::computeIndex(x, y);
131 }
132
155 IndexRange getIndices() const noexcept {
156 return IndexRange(typename IndexRange::iterator(), IndexRange::iterator::makeEnd(getOrder()));
157 }
158
161
181 template <typename Vector>
182 double sumWith(geom::Point2D const & point, Vector const & coefficients,
183 Workspace & workspace, SumMode mode=SumMode::FAST) const {
184 assert(workspace.getOrder() >= getOrder());
185 _basis1d.fill(point.getX(), workspace._x);
186 _basis1d.fill(point.getY(), workspace._y);
187 // This universal lambda lets us effectively template most of the
188 // implementation of this function on double vs. SafeSum<double>
189 // without having to define an external template.
190 auto accumulate = [coefficients, &workspace, this](auto & sum) {
191 for (auto const & index : getIndices()) {
192 sum += coefficients[index.flat]*workspace._x[index.nx]*workspace._y[index.ny];
193 }
194 };
195 double result = 0.0;
196 if (mode == SumMode::FAST) {
197 double z = 0.0;
198 accumulate(z);
199 result = z;
200 } else {
202 accumulate(z);
203 result = static_cast<double>(z);
204 }
205 return result;
206 }
207
209 template <typename Vector>
210 double sumWith(geom::Point2D const & point, Vector const & coefficients,
211 SumMode mode=SumMode::FAST) const {
212 auto workspace = makeWorkspace();
213 return sumWith(point, coefficients, workspace, mode);
214 }
215
227 template <typename Vector>
228 void fill(geom::Point2D const & point, Vector && basis, Workspace & workspace) const {
229 assert(workspace.getOrder() >= getOrder());
230 _basis1d.fill(point.getX(), workspace._x);
231 _basis1d.fill(point.getY(), workspace._y);
232 for (auto const & index : getIndices()) {
233 std::forward<Vector>(basis)[index.flat] = workspace._x[index.nx]*workspace._y[index.ny];
234 }
235 }
236
238 template <typename Vector>
239 void fill(geom::Point2D const & point, Vector && basis) const {
240 auto workspace = makeWorkspace();
241 fill(point, std::forward<Vector>(basis), workspace);
242 }
243
244private:
245 Basis1d _basis1d;
246};
247
248}}} // namespace lsst::geom::polynomials
249
250#endif // !LSST_AFW_MATH_POLYNOMIALS_PackedBasis2d_h_INCLUDED
A basis interface for 1-d series expansions.
Definition Basis1d.h:36
A 2-d function defined by a series expansion and its coefficients.
Definition Function2d.h:42
A Basis2d formed from the product of a Basis1d for each of x and y, truncated at the sum of their ord...
void fill(geom::Point2D const &point, Vector &&basis, Workspace &workspace) const
double sumWith(geom::Point2D const &point, Vector const &coefficients, SumMode mode=SumMode::FAST) const
void fill(geom::Point2D const &point, Vector &&basis) const
static constexpr std::size_t computeSize(std::size_t order)
PackedBasis2d(Basis1d const &basis1d)
Construct from a 1-d basis that will be used for both x and y.
std::size_t index(std::size_t x, std::size_t y) const
double sumWith(geom::Point2D const &point, Vector const &coefficients, Workspace &workspace, SumMode mode=SumMode::FAST) const
PackedBasis2d & operator=(PackedBasis2d const &)=default
A workspace object that can be used to avoid extra memory allocations in repeated calls to PackedBasi...
std::size_t getOrder() const
Return the maximum order this workspace supports.
PackedBasisWorkspace2d(std::size_t order)
Construct workspace for a basis with the given order.
A specialized iterator range class for PackedIndexIterator, providing size calculation,...
PackedIndexIterator< packing > iterator
static constexpr std::size_t computeSize(std::size_t order) noexcept
Return the flattened size of an expansion with the given maximum order (inclusive).
static constexpr std::size_t computeIndex(std::size_t nx, std::size_t ny) noexcept
Return the flattened index for the element with the given x and y orders.
A numerically stable summation algorithm for floating-point numbers.
Definition SafeSum.h:62
A 2-d basis that transforms all input points before evaluating nested basis.
A 2-d separable affine transform that can be used to map one interval to another.
Definition Scaling2d.h:48
T forward(T... args)
Low-level polynomials (including special polynomials) in C++.
SumMode
Enum used to control how to sum polynomial terms.
Definition SafeSum.h:32
@ FAST
Summation using regular floating-point addition.
Definition SafeSum.h:34
Point< double, 2 > Point2D
Definition Point.h:324
STL namespace.