LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
aggregates.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2/*
3 * LSST Data Management System
4 * Copyright 2008-2014 LSST Corporation.
5 *
6 * This product includes software developed by the
7 * LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20 * the GNU General Public License along with this program. If not,
21 * see <http://www.lsstcorp.org/LegalNotices/>.
22 */
23
25#include "lsst/geom/Box.h"
28
29namespace lsst {
30namespace afw {
31namespace table {
32
33//============ PointKey =====================================================================================
34
35template <typename T>
37 std::string const &unit) {
38 Key<T> xKey = schema.addField<T>(schema.join(name, "x"), doc, unit);
39 Key<T> yKey = schema.addField<T>(schema.join(name, "y"), doc, unit);
40 return PointKey<T>(xKey, yKey);
41}
42
43template <typename T>
45 return lsst::geom::Point<T, 2>(record.get(_x), record.get(_y));
46}
47
48template <typename T>
49void PointKey<T>::set(BaseRecord &record, lsst::geom::Point<T, 2> const &value) const {
50 record.set(_x, value.getX());
51 record.set(_y, value.getY());
52}
53
54template class PointKey<int>;
55template class PointKey<double>;
56
57//============ Point3Key ====================================================================================
58
59template <typename T>
61 std::string const &unit) {
62 Key<T> xKey = schema.addField<T>(schema.join(name, "x"), doc, unit);
63 Key<T> yKey = schema.addField<T>(schema.join(name, "y"), doc, unit);
64 Key<T> zKey = schema.addField<T>(schema.join(name, "z"), doc, unit);
65 return Point3Key<T>(xKey, yKey, zKey);
66}
67
68template <typename T>
70 return lsst::geom::Point<T, 3>(record.get(_x), record.get(_y), record.get(_z));
71}
72
73template <typename T>
74void Point3Key<T>::set(BaseRecord &record, lsst::geom::Point<T, 3> const &value) const {
75 record.set(_x, value.getX());
76 record.set(_y, value.getY());
77 record.set(_z, value.getZ());
78}
79
80template class Point3Key<int>;
81template class Point3Key<double>;
82
83//============ BoxKey =====================================================================================
84
85template <typename Box>
87 std::string const &unit) {
88 auto minKey = PointKey<Element>::addFields(schema, schema.join(name, "min"), doc + " (minimum)", unit);
89 auto maxKey = PointKey<Element>::addFields(schema, schema.join(name, "max"), doc + " (maximum)", unit);
90 return BoxKey<Box>(minKey, maxKey);
92
93template <typename Box>
94Box BoxKey<Box>::get(BaseRecord const &record) const {
95 return Box(record.get(_min), record.get(_max), /*invert=*/false);
96}
97
98template <typename Box>
99void BoxKey<Box>::set(BaseRecord &record, Box const &value) const {
100 _min.set(record, value.getMin());
101 _max.set(record, value.getMax());
102}
103
104template class BoxKey<lsst::geom::Box2I>;
105template class BoxKey<lsst::geom::Box2D>;
106
107//============ CoordKey =====================================================================================
108
109CoordKey CoordKey::addFields(Schema &schema, std::string const &name, std::string const &doc) {
110 Key<lsst::geom::Angle> ra = schema.addField<lsst::geom::Angle>(schema.join(name, "ra"), doc);
111 Key<lsst::geom::Angle> dec = schema.addField<lsst::geom::Angle>(schema.join(name, "dec"), doc);
112 return CoordKey(ra, dec);
113}
114
116 return lsst::geom::SpherePoint(record.get(_ra), record.get(_dec));
117}
118
119void CoordKey::set(BaseRecord &record, lsst::geom::SpherePoint const &value) const {
120 record.set(_ra, value.getLongitude());
121 record.set(_dec, value.getLatitude());
122}
123
125 return ErrorKey(schema["coord"], {"ra", "dec"});
126}
127
129 return ErrorKey::addFields(schema, "coord", {"ra", "dec"}, "rad");
130}
131
132//============ QuadrupoleKey ================================================================================
133
136 std::string unit = coordType == CoordinateType::PIXEL ? "pixel^2" : "rad^2";
137
138 Key<double> xxKey = schema.addField<double>(schema.join(name, "xx"), doc, unit);
139 Key<double> yyKey = schema.addField<double>(schema.join(name, "yy"), doc, unit);
140 Key<double> xyKey = schema.addField<double>(schema.join(name, "xy"), doc, unit);
141 return QuadrupoleKey(xxKey, yyKey, xyKey);
142}
143
145 return geom::ellipses::Quadrupole(record.get(_ixx), record.get(_iyy), record.get(_ixy));
146}
147
149 record.set(_ixx, value.getIxx());
150 record.set(_iyy, value.getIyy());
151 record.set(_ixy, value.getIxy());
152}
153
154//============ EllipseKey ================================================================================
155
157 std::string const &unit) {
159 PointKey<double> pKey = PointKey<double>::addFields(schema, name, doc, unit);
160 return EllipseKey(qKey, pKey);
161}
162
164 return geom::ellipses::Ellipse(record.get(_qKey), record.get(_pKey));
165}
166
167void EllipseKey::set(BaseRecord &record, geom::ellipses::Ellipse const &value) const {
168 _qKey.set(record, value.getCore());
169 _pKey.set(record, value.getCenter());
170}
171
172//============ CovarianceMatrixKey ==========================================================================
173
174template <typename T, int N>
176 NameArray const &names,
177 std::string const &unit, bool diagonalOnly) {
178 NameArray units(names.size(), unit);
179 return addFields(schema, prefix, names, units, diagonalOnly);
180}
181
182template <typename T, int N>
184 NameArray const &names, NameArray const &units,
185 bool diagonalOnly) {
186 if (N != Eigen::Dynamic) {
188 "Size of names array (%d) does not match template argument (%d)");
190 "Size of units array (%d) does not match template argument (%d)");
191 }
194 err.reserve(names.size());
195 for (std::size_t i = 0; i < names.size(); ++i) {
196 err.push_back(schema.addField<T>(schema.join(prefix, names[i] + "Err"),
197 "1-sigma uncertainty on " + names[i], units[i]));
198 }
199 if (!diagonalOnly) {
200 cov.reserve((names.size() * (names.size() - 1)) / 2);
201 for (std::size_t i = 0; i < names.size(); ++i) {
202 for (std::size_t j = 0; j < i; ++j) {
203 // We iterate over the lower-triangular part of the matrix in row-major order,
204 // but we use the upper-triangular names (i.e. we switch the order of i and j, below).
205 // That puts the elements in the order expected by the constructor we call below,
206 // while creating the field names users would expect from the ordering of their name
207 // vector (i.e. _a_b_Cov instead of _b_a_Cov if names=[a, b]).
208 cov.push_back(schema.addField<T>(
209 schema.join(prefix, names[j], names[i], "Cov"),
210 "uncertainty covariance between " + names[j] + " and " + names[i],
211 units[j] + (units[j].empty() || units[i].empty() ? "" : " ") + units[i]));
212 }
213 }
214 }
216}
217
218template <typename T, int N>
220
221template <typename T, int N>
223 : _err(err), _cov(cov) {
224 if (N != Eigen::Dynamic) {
226 "Size of err array (%d) does not match template argument (%d)");
227 }
228 if (!cov.empty()) {
229 LSST_THROW_IF_NE(cov.size(), err.size() * (err.size() - 1) / 2, pex::exceptions::LengthError,
230 "Size of cov array (%d) is does not match with size inferred from err array (%d)");
231 bool haveCov = false;
232 for (typename CovarianceKeyArray::const_iterator i = _cov.begin(); i != _cov.end(); ++i) {
233 if (i->isValid()) haveCov = true;
234 }
235 if (!haveCov) _cov.resize(0);
236 }
237}
238
239template <typename T, int N>
241 : _err(names.size()), _cov(names.size() * (names.size() - 1) / 2) {
242 int const n = names.size();
243 int k = 0;
244 bool haveCov = false;
245 for (int i = 0; i < n; ++i) {
246 _err[i] = s[names[i] + "Err"];
247 for (int j = 0; j < i; ++j, ++k) {
248 try {
249 _cov[k] = s[names[i] + "_" + names[j] + "_Cov"];
250 haveCov = true;
252 try {
253 _cov[k] = s[names[j] + "_" + names[i] + "_Cov"];
254 haveCov = true;
256 }
257 }
259 }
260 if (!haveCov) _cov.resize(0);
262
263template <typename T, int N>
265template <typename T, int N>
267template <typename T, int N>
269template <typename T, int N>
271template <typename T, int N>
273
274// these are workarounds for the fact that Eigen has different constructors for
275// dynamic-sized matrices and fixed-size matrices, but we don't want to have to
276// partial-specialize the entire template just to change one line
277namespace {
278
279template <typename T, int N>
280Eigen::Matrix<T, N, N> makeZeroMatrix(int n, CovarianceMatrixKey<T, N> const *) {
281 return Eigen::Matrix<T, N, N>::Zero();
282}
283
284template <typename T>
285Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> makeZeroMatrix(
286 int n, CovarianceMatrixKey<T, Eigen::Dynamic> const *) {
287 return Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(n, n);
288}
289
290} // namespace
291
292template <typename T, int N>
293Eigen::Matrix<T, N, N> CovarianceMatrixKey<T, N>::get(BaseRecord const &record) const {
294 Eigen::Matrix<T, N, N> value = makeZeroMatrix(_err.size(), this);
295 int const n = _err.size();
296 int k = 0;
297 for (int i = 0; i < n; ++i) {
298 T err = record.get(_err[i]);
299 value(i, i) = err * err;
300 if (!_cov.empty()) {
301 for (int j = 0; j < i; ++j, ++k) {
302 if (_cov[k].isValid()) {
303 value(i, j) = value(j, i) = record.get(_cov[k]);
304 }
305 }
306 }
307 }
308 return value;
309}
310
311template <typename T, int N>
312void CovarianceMatrixKey<T, N>::set(BaseRecord &record, Eigen::Matrix<T, N, N> const &value) const {
313 int const n = _err.size();
314 int k = 0;
315 for (int i = 0; i < n; ++i) {
316 record.set(_err[i], std::sqrt(value(i, i)));
317 if (!_cov.empty()) {
318 for (int j = 0; j < i; ++j, ++k) {
319 if (_cov[k].isValid()) {
320 record.set(_cov[k], value(i, j));
321 }
322 }
323 }
324 }
325}
326
327template <typename T, int N>
329 int const n = _err.size();
330 if (n < 1) return false;
331 for (int i = 0; i < n; ++i) {
332 if (!_err[i].isValid()) return false;
333 }
334 return true;
335}
336
337template <typename T, int N>
339 if (_err.size() != other._err.size()) {
340 return false;
341 }
342 if (_cov.size() != other._cov.size()) {
343 return false;
344 }
345 int const n = _err.size();
346 int k = 0;
347 for (int i = 0; i < n; ++i) {
348 if (_err[i] != other._err[i]) {
349 return false;
350 }
351 if (!_cov.empty()) {
352 for (int j = 0; j < i; ++j, ++k) {
353 if (_cov[k] != other._cov[k]) {
354 return false;
355 }
356 }
357 }
358 }
359 return true;
360}
361
362template <typename T, int N>
364 // Completely arbitrary seeds, different to avoid any weird degeneracies/interactions
365 return utils::hashCombine(17, utils::hashIterable(19, _err.begin(), _err.end()),
366 utils::hashIterable(23, _cov.begin(), _cov.end()));
367}
368
369template <typename T, int N>
370T CovarianceMatrixKey<T, N>::getElement(BaseRecord const &record, int i, int j) const {
371 if (i == j) {
372 T err = record.get(_err[i]);
373 return err * err;
374 }
375 if (_cov.empty()) {
376 return 0.0;
377 }
378 Key<T> key = (i < j) ? _cov[j * (j - 1) / 2 + i] : _cov[i * (i - 1) / 2 + j];
379 return key.isValid() ? record.get(key) : 0.0;
380}
381
382template <typename T, int N>
383void CovarianceMatrixKey<T, N>::setElement(BaseRecord &record, int i, int j, T value) const {
384 if (i == j) {
385 record.set(_err[i], std::sqrt(value));
386 } else {
387 if (_cov.empty()) {
388 throw LSST_EXCEPT(
390 (boost::format("Cannot set covariance element %d,%d; no fields for covariance") % i % j)
391 .str());
392 }
393 Key<T> key = (i < j) ? _cov[j * (j - 1) / 2 + i] : _cov[i * (i - 1) / 2 + j];
394 if (!key.isValid()) {
395 throw LSST_EXCEPT(
397 (boost::format("Cannot set covariance element %d,%d; no field for this element") % i % j)
398 .str());
399 }
400 record.set(key, value);
401 }
402}
403
404template class CovarianceMatrixKey<float, 2>;
405template class CovarianceMatrixKey<float, 3>;
406template class CovarianceMatrixKey<float, 4>;
407template class CovarianceMatrixKey<float, 5>;
409template class CovarianceMatrixKey<double, 2>;
410template class CovarianceMatrixKey<double, 3>;
411template class CovarianceMatrixKey<double, 4>;
412template class CovarianceMatrixKey<double, 5>;
414} // namespace table
415} // namespace afw
416} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
double dec
Definition Match.cc:41
std::string prefix
#define LSST_THROW_IF_NE(N1, N2, EXC_CLASS, MSG)
Check whether the given values are equal, and throw an LSST Exception if they are not.
Definition asserts.h:38
T begin(T... args)
An ellipse defined by an arbitrary BaseCore and a center point.
Definition Ellipse.h:51
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
Tag types used to declare specialized field types.
Definition misc.h:31
Base class for all records.
Definition BaseRecord.h:31
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition BaseRecord.h:151
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition BaseRecord.h:164
Box get(BaseRecord const &record) const override
Get a Box from the given record.
Definition aggregates.cc:94
static BoxKey addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add _min_x, _min_y, _max_x, _max_y fields to a Schema, and return a BoxKey that points to them.
Definition aggregates.cc:86
void set(BaseRecord &record, Box const &value) const override
Set a Box in the given record.
Definition aggregates.cc:99
A FunctorKey used to get or set celestial coordinates from a pair of lsst::geom::Angle keys.
Definition aggregates.h:292
CovarianceMatrixKey< float, 2 > ErrorKey
Definition aggregates.h:304
static CoordKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc)
Add a pair of _ra, _dec fields to a Schema, and return a CoordKey that points to them.
CoordKey() noexcept
Default constructor; instance will not be usable unless subsequently assigned to.
Definition aggregates.h:310
void set(BaseRecord &record, lsst::geom::SpherePoint const &value) const override
Set an lsst::geom::SpherePoint in the given record.
static ErrorKey addErrorFields(Schema &schema)
static ErrorKey getErrorKey(Schema const &schema)
lsst::geom::SpherePoint get(BaseRecord const &record) const override
Get an lsst::geom::SpherePoint from the given record.
bool operator==(CovarianceMatrixKey const &other) const noexcept
Compare the FunctorKey for equality with another, using its constituent Keys.
CovarianceMatrixKey()
Construct an invalid instance; must assign before subsequent use.
Eigen::Matrix< T, N, N > get(BaseRecord const &record) const override
Get a covariance matrix from the given record.
void setElement(BaseRecord &record, int i, int j, T value) const
Set the element in row i and column j.
T getElement(BaseRecord const &record, int i, int j) const
Return the element in row i and column j.
void set(BaseRecord &record, Eigen::Matrix< T, N, N > const &value) const override
Set a covariance matrix in the given record (uses only the lower triangle of the given matrix)
~CovarianceMatrixKey() noexcept override
static CovarianceMatrixKey addFields(Schema &schema, std::string const &prefix, NameArray const &names, std::string const &unit, bool diagonalOnly=false)
Add covariance matrix fields to a Schema, and return a CovarianceMatrixKey to manage them.
std::size_t hash_value() const noexcept
Return a hash of this object.
bool isValid() const noexcept
Return True if all the constituent error Keys are valid.
CovarianceMatrixKey & operator=(CovarianceMatrixKey const &)
A FunctorKey used to get or set a geom::ellipses::Ellipse from an (xx,yy,xy,x,y) tuple of Keys.
Definition aggregates.h:447
void set(BaseRecord &record, geom::ellipses::Ellipse const &value) const override
Set an Ellipse in the given record.
static EllipseKey addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add a set of _xx, _yy, _xy, _x, _y fields to a Schema, and return an EllipseKey that points to them.
geom::ellipses::Ellipse get(BaseRecord const &record) const override
Get an Ellipse from the given record.
EllipseKey() noexcept
Default constructor; instance will not be usable unless subsequently assigned to.
Definition aggregates.h:463
void set(BaseRecord &record, lsst::geom::Point< T, 3 > const &value) const override
Set a Point in the given record.
Definition aggregates.cc:74
static Point3Key addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add a pair of _x, _y, _z fields to a Schema, and return a Point3Key that points to them.
Definition aggregates.cc:60
lsst::geom::Point< T, 3 > get(BaseRecord const &record) const override
Get a Point from the given record.
Definition aggregates.cc:69
void set(BaseRecord &record, lsst::geom::Point< T, 2 > const &value) const override
Set a Point in the given record.
Definition aggregates.cc:49
static PointKey addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add a pair of _x, _y fields to a Schema, and return a PointKey that points to them.
Definition aggregates.cc:36
lsst::geom::Point< T, 2 > get(BaseRecord const &record) const override
Get a Point from the given record.
Definition aggregates.cc:44
A FunctorKey used to get or set a geom::ellipses::Quadrupole from a tuple of constituent Keys.
Definition aggregates.h:369
geom::ellipses::Quadrupole get(BaseRecord const &record) const override
Get a Quadrupole from the given record.
QuadrupoleKey() noexcept
Default constructor; instance will not be usable unless subsequently assigned to.
Definition aggregates.h:386
void set(BaseRecord &record, geom::ellipses::Quadrupole const &value) const override
Set a Quadrupole in the given record.
static QuadrupoleKey addFields(Schema &schema, std::string const &name, std::string const &doc, CoordinateType coordType=CoordinateType::PIXEL)
Add a set of quadrupole subfields to a schema and return a QuadrupoleKey that points to them.
Defines the fields and offsets for a table.
Definition Schema.h:51
A proxy type for name lookups in a Schema.
Definition Schema.h:367
A class representing an angle.
Definition Angle.h:128
A coordinate class intended to represent absolute positions.
Definition Point.h:169
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Reports attempts to access elements using an invalid key.
Definition Runtime.h:151
T end(T... args)
bool isValid
Definition fits.cc:404
CoordinateType
Enum used to set units for geometric FunctorKeys.
Definition aggregates.h:364
std::size_t hashCombine(std::size_t seed) noexcept
Combine hashes.
Definition hashCombine.h:35
std::size_t hashIterable(std::size_t seed, InputIterator begin, InputIterator end) noexcept
Combine hashes in an iterable.
Definition hashCombine.h:93
T resize(T... args)
T sqrt(T... args)