LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
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 
29 namespace lsst {
30 namespace afw {
31 namespace table {
32 
33 //============ PointKey =====================================================================================
34 
35 template <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 
43 template <typename T>
45  return lsst::geom::Point<T, 2>(record.get(_x), record.get(_y));
46 }
47 
48 template <typename T>
49 void 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 
54 template class PointKey<int>;
55 template class PointKey<double>;
56 
57 //============ BoxKey =====================================================================================
58 
59 template <typename Box>
61  std::string const &unit) {
62  auto minKey = PointKey<Element>::addFields(schema, schema.join(name, "min"), doc + " (minimum)", unit);
63  auto maxKey = PointKey<Element>::addFields(schema, schema.join(name, "max"), doc + " (maximum)", unit);
64  return BoxKey<Box>(minKey, maxKey);
65 }
66 
67 template <typename Box>
68 Box BoxKey<Box>::get(BaseRecord const &record) const {
69  return Box(record.get(_min), record.get(_max), /*invert=*/false);
70 }
71 
72 template <typename Box>
73 void BoxKey<Box>::set(BaseRecord &record, Box const &value) const {
74  _min.set(record, value.getMin());
75  _max.set(record, value.getMax());
76 }
77 
78 template class BoxKey<lsst::geom::Box2I>;
79 template class BoxKey<lsst::geom::Box2D>;
80 
81 //============ CoordKey =====================================================================================
82 
84  Key<lsst::geom::Angle> ra = schema.addField<lsst::geom::Angle>(schema.join(name, "ra"), doc);
85  Key<lsst::geom::Angle> dec = schema.addField<lsst::geom::Angle>(schema.join(name, "dec"), doc);
86  return CoordKey(ra, dec);
87 }
88 
90  return lsst::geom::SpherePoint(record.get(_ra), record.get(_dec));
91 }
92 
93 void CoordKey::set(BaseRecord &record, lsst::geom::SpherePoint const &value) const {
94  record.set(_ra, value.getLongitude());
95  record.set(_dec, value.getLatitude());
96 }
97 
98 //============ QuadrupoleKey ================================================================================
99 
101  CoordinateType coordType) {
102  std::string unit = coordType == CoordinateType::PIXEL ? "pixel^2" : "rad^2";
103 
104  Key<double> xxKey = schema.addField<double>(schema.join(name, "xx"), doc, unit);
105  Key<double> yyKey = schema.addField<double>(schema.join(name, "yy"), doc, unit);
106  Key<double> xyKey = schema.addField<double>(schema.join(name, "xy"), doc, unit);
107  return QuadrupoleKey(xxKey, yyKey, xyKey);
108 }
109 
111  return geom::ellipses::Quadrupole(record.get(_ixx), record.get(_iyy), record.get(_ixy));
112 }
113 
114 void QuadrupoleKey::set(BaseRecord &record, geom::ellipses::Quadrupole const &value) const {
115  record.set(_ixx, value.getIxx());
116  record.set(_iyy, value.getIyy());
117  record.set(_ixy, value.getIxy());
118 }
119 
120 //============ EllipseKey ================================================================================
121 
123  std::string const &unit) {
126  return EllipseKey(qKey, pKey);
127 }
128 
130  return geom::ellipses::Ellipse(record.get(_qKey), record.get(_pKey));
131 }
132 
133 void EllipseKey::set(BaseRecord &record, geom::ellipses::Ellipse const &value) const {
134  _qKey.set(record, value.getCore());
135  _pKey.set(record, value.getCenter());
136 }
137 
138 //============ CovarianceMatrixKey ==========================================================================
139 
140 template <typename T, int N>
142  NameArray const &names,
143  std::string const &unit, bool diagonalOnly) {
144  NameArray units(names.size(), unit);
145  return addFields(schema, prefix, names, units, diagonalOnly);
146 }
147 
148 template <typename T, int N>
150  NameArray const &names, NameArray const &units,
151  bool diagonalOnly) {
152  if (N != Eigen::Dynamic) {
154  "Size of names array (%d) does not match template argument (%d)");
156  "Size of units array (%d) does not match template argument (%d)");
157  }
158  ErrKeyArray err;
159  CovarianceKeyArray cov;
160  err.reserve(names.size());
161  for (std::size_t i = 0; i < names.size(); ++i) {
162  err.push_back(schema.addField<T>(schema.join(prefix, names[i] + "Err"),
163  "1-sigma uncertainty on " + names[i], units[i]));
164  }
165  if (!diagonalOnly) {
166  cov.reserve((names.size() * (names.size() - 1)) / 2);
167  for (std::size_t i = 0; i < names.size(); ++i) {
168  for (std::size_t j = 0; j < i; ++j) {
169  // We iterate over the lower-triangular part of the matrix in row-major order,
170  // but we use the upper-triangular names (i.e. we switch the order of i and j, below).
171  // That puts the elements in the order expected by the constructor we call below,
172  // while creating the field names users would expect from the ordering of their name
173  // vector (i.e. _a_b_Cov instead of _b_a_Cov if names=[a, b]).
174  cov.push_back(schema.addField<T>(
175  schema.join(prefix, names[j], names[i], "Cov"),
176  "uncertainty covariance between " + names[j] + " and " + names[i],
177  units[j] + (units[j].empty() || units[i].empty() ? "" : " ") + units[i]));
178  }
179  }
180  }
181  return CovarianceMatrixKey<T, N>(err, cov);
182 }
183 
184 template <typename T, int N>
186 
187 template <typename T, int N>
189  : _err(err), _cov(cov) {
190  if (N != Eigen::Dynamic) {
192  "Size of err array (%d) does not match template argument (%d)");
193  }
194  if (!cov.empty()) {
195  LSST_THROW_IF_NE(cov.size(), err.size() * (err.size() - 1) / 2, pex::exceptions::LengthError,
196  "Size of cov array (%d) is does not match with size inferred from err array (%d)");
197  bool haveCov = false;
198  for (typename CovarianceKeyArray::const_iterator i = _cov.begin(); i != _cov.end(); ++i) {
199  if (i->isValid()) haveCov = true;
200  }
201  if (!haveCov) _cov.resize(0);
202  }
203 }
204 
205 template <typename T, int N>
207  : _err(names.size()), _cov(names.size() * (names.size() - 1) / 2) {
208  int const n = names.size();
209  int k = 0;
210  bool haveCov = false;
211  for (int i = 0; i < n; ++i) {
212  _err[i] = s[names[i] + "Err"];
213  for (int j = 0; j < i; ++j, ++k) {
214  try {
215  _cov[k] = s[names[i] + "_" + names[j] + "_Cov"];
216  haveCov = true;
217  } catch (pex::exceptions::NotFoundError &) {
218  try {
219  _cov[k] = s[names[j] + "_" + names[i] + "_Cov"];
220  haveCov = true;
221  } catch (pex::exceptions::NotFoundError &) {
222  }
223  }
224  }
225  }
226  if (!haveCov) _cov.resize(0);
227 }
228 
229 template <typename T, int N>
231 template <typename T, int N>
233 template <typename T, int N>
235 template <typename T, int N>
237 template <typename T, int N>
239 
240 // these are workarounds for the fact that Eigen has different constructors for
241 // dynamic-sized matrices and fixed-size matrices, but we don't want to have to
242 // partial-specialize the entire template just to change one line
243 namespace {
244 
245 template <typename T, int N>
246 Eigen::Matrix<T, N, N> makeZeroMatrix(int n, CovarianceMatrixKey<T, N> const *) {
247  return Eigen::Matrix<T, N, N>::Zero();
248 }
249 
250 template <typename T>
251 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> makeZeroMatrix(
252  int n, CovarianceMatrixKey<T, Eigen::Dynamic> const *) {
253  return Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(n, n);
254 }
255 
256 } // namespace
257 
258 template <typename T, int N>
259 Eigen::Matrix<T, N, N> CovarianceMatrixKey<T, N>::get(BaseRecord const &record) const {
260  Eigen::Matrix<T, N, N> value = makeZeroMatrix(_err.size(), this);
261  int const n = _err.size();
262  int k = 0;
263  for (int i = 0; i < n; ++i) {
264  T err = record.get(_err[i]);
265  value(i, i) = err * err;
266  if (!_cov.empty()) {
267  for (int j = 0; j < i; ++j, ++k) {
268  if (_cov[k].isValid()) {
269  value(i, j) = value(j, i) = record.get(_cov[k]);
270  }
271  }
272  }
273  }
274  return value;
275 }
276 
277 template <typename T, int N>
278 void CovarianceMatrixKey<T, N>::set(BaseRecord &record, Eigen::Matrix<T, N, N> const &value) const {
279  int const n = _err.size();
280  int k = 0;
281  for (int i = 0; i < n; ++i) {
282  record.set(_err[i], std::sqrt(value(i, i)));
283  if (!_cov.empty()) {
284  for (int j = 0; j < i; ++j, ++k) {
285  if (_cov[k].isValid()) {
286  record.set(_cov[k], value(i, j));
287  }
288  }
289  }
290  }
291 }
292 
293 template <typename T, int N>
294 bool CovarianceMatrixKey<T, N>::isValid() const noexcept {
295  int const n = _err.size();
296  if (n < 1) return false;
297  for (int i = 0; i < n; ++i) {
298  if (!_err[i].isValid()) return false;
299  }
300  return true;
301 }
302 
303 template <typename T, int N>
305  if (_err.size() != other._err.size()) {
306  return false;
307  }
308  if (_cov.size() != other._cov.size()) {
309  return false;
310  }
311  int const n = _err.size();
312  int k = 0;
313  for (int i = 0; i < n; ++i) {
314  if (_err[i] != other._err[i]) {
315  return false;
316  }
317  if (!_cov.empty()) {
318  for (int j = 0; j < i; ++j, ++k) {
319  if (_cov[k] != other._cov[k]) {
320  return false;
321  }
322  }
323  }
324  }
325  return true;
326 }
327 
328 template <typename T, int N>
330  // Completely arbitrary seeds, different to avoid any weird degeneracies/interactions
331  return utils::hashCombine(17, utils::hashIterable(19, _err.begin(), _err.end()),
332  utils::hashIterable(23, _cov.begin(), _cov.end()));
333 }
334 
335 template <typename T, int N>
336 T CovarianceMatrixKey<T, N>::getElement(BaseRecord const &record, int i, int j) const {
337  if (i == j) {
338  T err = record.get(_err[i]);
339  return err * err;
340  }
341  if (_cov.empty()) {
342  return 0.0;
343  }
344  Key<T> key = (i < j) ? _cov[j * (j - 1) / 2 + i] : _cov[i * (i - 1) / 2 + j];
345  return key.isValid() ? record.get(key) : 0.0;
346 }
347 
348 template <typename T, int N>
349 void CovarianceMatrixKey<T, N>::setElement(BaseRecord &record, int i, int j, T value) const {
350  if (i == j) {
351  record.set(_err[i], std::sqrt(value));
352  } else {
353  if (_cov.empty()) {
354  throw LSST_EXCEPT(
356  (boost::format("Cannot set covariance element %d,%d; no fields for covariance") % i % j)
357  .str());
358  }
359  Key<T> key = (i < j) ? _cov[j * (j - 1) / 2 + i] : _cov[i * (i - 1) / 2 + j];
360  if (!key.isValid()) {
361  throw LSST_EXCEPT(
363  (boost::format("Cannot set covariance element %d,%d; no field for this element") % i % j)
364  .str());
365  }
366  record.set(key, value);
367  }
368 }
369 
370 template class CovarianceMatrixKey<float, 2>;
371 template class CovarianceMatrixKey<float, 3>;
372 template class CovarianceMatrixKey<float, 4>;
373 template class CovarianceMatrixKey<float, 5>;
375 template class CovarianceMatrixKey<double, 2>;
376 template class CovarianceMatrixKey<double, 3>;
377 template class CovarianceMatrixKey<double, 4>;
378 template class CovarianceMatrixKey<double, 5>;
380 } // namespace table
381 } // namespace afw
382 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
double dec
Definition: Match.cc:41
ItemVariant const * other
Definition: Schema.cc:56
Key< U > key
Definition: Schema.cc:281
std::string prefix
Definition: SchemaMapper.cc:79
table::Schema schema
Definition: python.h:134
#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
lsst::geom::Point2D const & getCenter() const
Return the center point.
Definition: Ellipse.h:62
BaseCore const & getCore() const
Return the ellipse core.
Definition: Ellipse.h:71
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:47
Base class for all records.
Definition: BaseRecord.h:31
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:164
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
A FunctorKey used to get or set a lsst::geom::Box2I or Box2D from a (min, max) pair of PointKeys.
Definition: aggregates.h:126
Box get(BaseRecord const &record) const override
Get a Box from the given record.
Definition: aggregates.cc:68
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:60
void set(BaseRecord &record, Box const &value) const override
Set a Box in the given record.
Definition: aggregates.cc:73
A FunctorKey used to get or set celestial coordinates from a pair of lsst::geom::Angle keys.
Definition: aggregates.h:210
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.
Definition: aggregates.cc:83
CoordKey() noexcept
Default constructor; instance will not be usable unless subsequently assigned to.
Definition: aggregates.h:223
void set(BaseRecord &record, lsst::geom::SpherePoint const &value) const override
Set an lsst::geom::SpherePoint in the given record.
Definition: aggregates.cc:93
lsst::geom::SpherePoint get(BaseRecord const &record) const override
Get an lsst::geom::SpherePoint from the given record.
Definition: aggregates.cc:89
bool operator==(CovarianceMatrixKey const &other) const noexcept
Compare the FunctorKey for equality with another, using its constituent Keys.
Definition: aggregates.cc:304
CovarianceMatrixKey()
Construct an invalid instance; must assign before subsequent use.
Definition: aggregates.cc:185
Eigen::Matrix< T, N, N > get(BaseRecord const &record) const override
Get a covariance matrix from the given record.
Definition: aggregates.cc:259
void setElement(BaseRecord &record, int i, int j, T value) const
Set the element in row i and column j.
Definition: aggregates.cc:349
T getElement(BaseRecord const &record, int i, int j) const
Return the element in row i and column j.
Definition: aggregates.cc:336
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)
Definition: aggregates.cc:278
~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.
Definition: aggregates.cc:141
std::size_t hash_value() const noexcept
Return a hash of this object.
Definition: aggregates.cc:329
bool isValid() const noexcept
Return True if all the constituent error Keys are valid.
Definition: aggregates.cc:294
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:360
void set(BaseRecord &record, geom::ellipses::Ellipse const &value) const override
Set an Ellipse in the given record.
Definition: aggregates.cc:133
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.
Definition: aggregates.cc:122
geom::ellipses::Ellipse get(BaseRecord const &record) const override
Get an Ellipse from the given record.
Definition: aggregates.cc:129
EllipseKey() noexcept
Default constructor; instance will not be usable unless subsequently assigned to.
Definition: aggregates.h:376
A class used as a handle to a particular field in a table.
Definition: Key.h:53
A FunctorKey used to get or set a lsst::geom::Point from an (x,y) pair of int or double Keys.
Definition: aggregates.h:49
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:282
geom::ellipses::Quadrupole get(BaseRecord const &record) const override
Get a Quadrupole from the given record.
Definition: aggregates.cc:110
QuadrupoleKey() noexcept
Default constructor; instance will not be usable unless subsequently assigned to.
Definition: aggregates.h:299
void set(BaseRecord &record, geom::ellipses::Quadrupole const &value) const override
Set a Quadrupole in the given record.
Definition: aggregates.cc:114
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.
Definition: aggregates.cc:100
Defines the fields and offsets for a table.
Definition: Schema.h:50
A proxy type for name lookups in a Schema.
Definition: Schema.h:357
A class representing an angle.
Definition: Angle.h:127
A coordinate class intended to represent absolute positions (2-d specialization).
Definition: Point.h:211
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
Angle getLatitude() const noexcept
The latitude of this point.
Definition: SpherePoint.h:190
Angle getLongitude() const noexcept
The longitude of this point.
Definition: SpherePoint.h:178
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 empty(T... args)
T end(T... args)
bool isValid
Definition: fits.cc:399
CoordinateType
Enum used to set units for geometric FunctorKeys.
Definition: aggregates.h:277
lsst::geom::SpherePoint SpherePoint
Definition: misc.h:35
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
std::size_t hashIterable(std::size_t seed, InputIterator begin, InputIterator end) noexcept
Combine hashes in an iterable.
Definition: hashCombine.h:93
std::size_t hashCombine(std::size_t seed) noexcept
Combine hashes.
Definition: hashCombine.h:35
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
T resize(T... args)
T size(T... args)
T sqrt(T... args)