LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
TransmissionCurve.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2017 LSST/AURA.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <http://www.lsstcorp.org/LegalNotices/>.
21  */
22 
23 #include <algorithm>
24 #include <memory>
25 
26 #include "ndarray.h"
27 
28 #include "gsl/gsl_interp.h"
29 #include "gsl/gsl_interp2d.h"
30 #include "gsl/gsl_errno.h"
31 
36 #include "lsst/afw/table/io/Persistable.cc"
37 
38 using namespace std::string_literals;
39 
40 #define LSST_CHECK_GSL(type, status) \
41  if (status) throw LSST_EXCEPT(type, "GSL error: "s + ::gsl_strerror(status))
42 
43 namespace lsst {
44 namespace afw {
45 
46 template std::shared_ptr<image::TransmissionCurve> table::io::PersistableFacade<
48 
49 namespace image {
50 
51 namespace {
52 
53 /*
54  * The TransmissionCurve implementation returned by TransmissionCurve::makeIdentity.
55  *
56  * This is zero-state singleton whose throughput is always exactly 1.
57  */
58 class IdentityTransmissionCurve : public TransmissionCurve {
59 public:
60  static constexpr char const* NAME = "IdentityTransmissionCurve";
61 
63  static std::shared_ptr<IdentityTransmissionCurve> instance(new IdentityTransmissionCurve());
64  return instance;
65  }
66 
67  std::pair<double, double> getWavelengthBounds() const override {
68  constexpr double inf = std::numeric_limits<double>::infinity();
69  return std::make_pair(-inf, inf);
70  }
71 
72  std::pair<double, double> getThroughputAtBounds() const override { return std::make_pair(1.0, 1.0); }
73 
74  void sampleAt(lsst::geom::Point2D const&, ndarray::Array<double const, 1, 1> const& wavelengths,
75  ndarray::Array<double, 1, 1> const& out) const override {
76  out.deep() = 1.0;
77  }
78 
79  bool isPersistable() const noexcept override { return true; }
80 
81 protected:
82  // transforming an IdentityTransmissionCurve is a no-op
85  return shared_from_this();
86  }
87 
88  // multiplying an IdentityTransmissionCurve always yields the other operand
91  return other;
92  }
93 
94  std::string getPersistenceName() const override { return NAME; }
95 
96  void write(OutputArchiveHandle& handle) const override { handle.saveEmpty(); }
97 
98  class Factory : public table::io::PersistableFactory {
99  public:
100  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
101  CatalogVector const& catalogs) const override {
102  LSST_ARCHIVE_ASSERT(catalogs.size() == 0u);
103  return get();
104  }
105 
106  Factory() : table::io::PersistableFactory(NAME) {}
107  };
108 
109  static Factory registration;
110 
111 private:
112  IdentityTransmissionCurve() = default;
113 };
114 
115 IdentityTransmissionCurve::Factory IdentityTransmissionCurve::registration;
116 
117 /*
118  * InterpolatedTransmissionCurve: implements makeSpatiallyConstant and makeRadial.
119  *
120  * InterpolatedTransmissionCurve is templated on an implementation class
121  * that does most of the work, letting it handle the boilerplate that is
122  * common to both spatially-constant and radially-varying curves.
123  *
124  * We use two tricks to avoid repetition that bear some explanation:
125  *
126  * - Even though the two implementations have different state, we can use
127  * the same PersistenceHelper template for both by using a std::vector
128  * to hold the keys for the arrays. Impl2d has one extra array to save,
129  * so its vector is one element longer.
130  *
131  * - The throughput array held by Impl2d is conceptually a 2-d array, but we
132  * actually store a flattened view into it. That's okay because both GSL
133  * and the persistence layer only care about the pointer and the size. That
134  * requires going through some hoops in makeRadial, but from then on it's
135  * clean sailing.
136  */
137 
138 template <typename T>
139 using GslPtr = std::unique_ptr<T, void (*)(T*)>;
140 
141 template <typename T>
142 GslPtr<T> makeGslPtr(T* p, void (*free)(T*)) {
143  if (p == nullptr) {
144  throw std::bad_alloc();
145  }
146  return GslPtr<T>(p, free);
147 }
148 
149 using ArrayKeyVector = std::vector<table::Key<table::Array<double>>>;
150 
151 /*
152  * Implementation object as a template parameter for the instantiation of
153  * InterpolatedTransmissionCurve used by makeSpatiallyConstant.
154  */
155 class Impl1d {
156 public:
157  static constexpr bool isSpatiallyConstant = true;
158  static constexpr char const* NAME = "SpatiallyConstantTransmissionCurve";
159 
160  // Initialize the GSL interpolator and take ownership of the given arrays.
161  // Array size consistency checks are done by makeSpatiallyConstant,
162  // so we don't need to worry about them here.
163  Impl1d(ndarray::Array<double, 1, 1> const& throughput, ndarray::Array<double, 1, 1> const& wavelengths)
164  : _throughput(throughput),
165  _wavelengths(wavelengths),
166  _interp(makeGslPtr(::gsl_interp_alloc(::gsl_interp_linear, _wavelengths.size()),
167  &::gsl_interp_free)) {
168  int status = ::gsl_interp_init(_interp.get(), _wavelengths.getData(), _throughput.getData(),
169  _wavelengths.size());
170  LSST_CHECK_GSL(pex::exceptions::LogicError, status);
171  }
172 
173  std::pair<double, double> getWavelengthBounds() const {
174  return std::make_pair(_wavelengths.front(), _wavelengths.back());
175  }
176 
177  static void setupPersistence(table::Schema& schema, ArrayKeyVector& keys) {
178  keys.push_back(schema.addField<table::Array<double>>("throughput", "array of known throughput values",
179  "", 0));
180  keys.push_back(schema.addField<table::Array<double>>(
181  "wavelengths", "array of known wavelength values", "angstrom", 0));
182  }
183 
184  void persist(table::BaseRecord& record, ArrayKeyVector const& keys) const {
185  record.set(keys[0], _throughput);
186  record.set(keys[1], _wavelengths);
187  }
188 
189  static Impl1d unpersist(table::BaseRecord& record, ArrayKeyVector const& keys) {
190  return Impl1d(record[keys[0]], record[keys[1]]);
191  }
192 
193  // A helper object constructed every time InterpolatedTransmissionCurve::sampleAt
194  // is called, and then invoked at every iteration of the loop therein.
195  struct Functor {
196  Functor(Impl1d const&, lsst::geom::Point2D const&)
197  : _accel(makeGslPtr(::gsl_interp_accel_alloc(), &gsl_interp_accel_free)) {}
198 
199  double operator()(Impl1d const& parent, double wavelength) {
200  double result = 0.0;
201  int status = ::gsl_interp_eval_e(parent._interp.get(), parent._wavelengths.getData(),
202  parent._throughput.getData(), wavelength, _accel.get(), &result);
203  LSST_CHECK_GSL(pex::exceptions::RuntimeError, status);
204  return result;
205  }
206 
207  private:
208  GslPtr<::gsl_interp_accel> _accel;
209  };
210 
211 private:
212  ndarray::Array<double, 1, 1> _throughput;
213  ndarray::Array<double, 1, 1> _wavelengths;
214  GslPtr<::gsl_interp> _interp;
215 };
216 
217 /*
218  * Implementation object as a template parameter for the instantiation of
219  * InterpolatedTransmissionCurve used by makeRadial.
220  */
221 class Impl2d {
222 public:
223  static constexpr bool isSpatiallyConstant = false;
224  static constexpr char const* NAME = "RadialTransmissionCurve";
225 
226  Impl2d(ndarray::Array<double, 1, 1> const& throughput, ndarray::Array<double, 1, 1> const& wavelengths,
227  ndarray::Array<double, 1, 1> const& radii)
228  : _throughput(throughput),
229  _wavelengths(wavelengths),
230  _radii(radii),
231  _interp(makeGslPtr(
232  ::gsl_interp2d_alloc(::gsl_interp2d_bilinear, _wavelengths.size(), _radii.size()),
233  &::gsl_interp2d_free)) {
234  int status = ::gsl_interp2d_init(_interp.get(), _wavelengths.getData(), _radii.getData(),
235  _throughput.getData(), _wavelengths.size(), _radii.size());
236  LSST_CHECK_GSL(pex::exceptions::LogicError, status);
237  }
238 
239  std::pair<double, double> getWavelengthBounds() const {
240  return std::make_pair(_wavelengths.front(), _wavelengths.back());
241  }
242 
243  static void setupPersistence(table::Schema& schema, ArrayKeyVector& keys) {
244  keys.push_back(schema.addField<table::Array<double>>(
245  "throughput",
246  "flattenned 2-d array of known throughput values, with radius dimension fastest", "", 0));
247  keys.push_back(schema.addField<table::Array<double>>(
248  "wavelengths", "array of known wavelength values", "angstrom", 0));
249  keys.push_back(schema.addField<table::Array<double>>("radii", "array of known radius values", "", 0));
250  }
251 
252  void persist(table::BaseRecord& record, ArrayKeyVector const& keys) const {
253  record.set(keys[0], _throughput);
254  record.set(keys[1], _wavelengths);
255  record.set(keys[2], _radii);
256  }
257 
258  static Impl2d unpersist(table::BaseRecord& record, ArrayKeyVector const& keys) {
259  return Impl2d(record[keys[0]], record[keys[1]], record[keys[2]]);
260  }
261 
262  // A helper object constructed every time InterpolatedTransmissionCurve::sampleAt
263  // is called, and then invoked at every iteration of the loop therein.
264  struct Functor {
265  Functor(Impl2d const& parent, lsst::geom::Point2D const& point)
266  : _radius(point.asEigen().norm()),
267  _radiusAccel(makeGslPtr(::gsl_interp_accel_alloc(), &gsl_interp_accel_free)),
268  _wavelengthAccel(makeGslPtr(::gsl_interp_accel_alloc(), &gsl_interp_accel_free)) {
269  _radius = std::max(_radius, parent._radii.front());
270  _radius = std::min(_radius, parent._radii.back());
271  }
272 
273  double operator()(Impl2d const& parent, double wavelength) {
274  double result = 0.0;
275  int status =
276  ::gsl_interp2d_eval_e(parent._interp.get(), parent._wavelengths.getData(),
277  parent._radii.getData(), parent._throughput.getData(), wavelength,
278  _radius, _radiusAccel.get(), _wavelengthAccel.get(), &result);
279  LSST_CHECK_GSL(pex::exceptions::RuntimeError, status);
280  return result;
281  }
282 
283  private:
284  double _radius;
285  GslPtr<::gsl_interp_accel> _radiusAccel;
286  GslPtr<::gsl_interp_accel> _wavelengthAccel;
287  };
288 
289 private:
290  ndarray::Array<double, 1, 1> _throughput;
291  ndarray::Array<double, 1, 1> _wavelengths;
292  ndarray::Array<double, 1, 1> _radii;
293  GslPtr<::gsl_interp2d> _interp;
294 };
295 
296 template <typename Impl>
297 class InterpolatedTransmissionCurve : public TransmissionCurve {
298 public:
299  InterpolatedTransmissionCurve(Impl impl, std::pair<double, double> throughputAtBounds)
300  : _atBounds(throughputAtBounds), _impl(std::move(impl)) {
301  if (!std::isfinite(_atBounds.first) || !std::isfinite(_atBounds.second)) {
302  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
303  "Throughput values at bounds must be finite");
304  }
305  }
306 
307  std::pair<double, double> getWavelengthBounds() const override { return _impl.getWavelengthBounds(); }
308 
309  std::pair<double, double> getThroughputAtBounds() const override { return _atBounds; }
310 
311  void sampleAt(lsst::geom::Point2D const& point, ndarray::Array<double const, 1, 1> const& wavelengths,
312  ndarray::Array<double, 1, 1> const& out) const override {
313  LSST_THROW_IF_NE(wavelengths.getSize<0>(), out.getSize<0>(), pex::exceptions::LengthError,
314  "Length of wavelength array (%d) does not match size of output array (%d)");
315  typename Impl::Functor functor(_impl, point);
316  auto bounds = _impl.getWavelengthBounds();
317  auto wlIter = wavelengths.begin();
318  for (auto outIter = out.begin(); outIter != out.end(); ++outIter, ++wlIter) {
319  double& y = *outIter;
320  if (*wlIter < bounds.first) {
321  y = _atBounds.first;
322  } else if (*wlIter > bounds.second) {
323  y = _atBounds.second;
324  } else {
325  y = functor(_impl, *wlIter);
326  }
327  }
328  }
329 
330  bool isPersistable() const noexcept override { return true; }
331 
332 protected:
333  std::shared_ptr<TransmissionCurve const> _transformedByImpl(
335  if (_impl.isSpatiallyConstant) {
336  return shared_from_this();
337  } else {
338  return TransmissionCurve::_transformedByImpl(transform);
339  }
340  }
341 
342  std::string getPersistenceName() const override { return Impl::NAME; }
343 
344  struct PersistenceHelper {
345  table::Schema schema;
346  table::Key<double> throughputAtMin;
347  table::Key<double> throughputAtMax;
348  ArrayKeyVector arrays;
349 
350  static PersistenceHelper const& get() {
351  static PersistenceHelper const instance;
352  return instance;
353  }
354 
355  private:
356  PersistenceHelper()
357  : schema(),
359  schema.addField<double>("throughputAtMin", "throughput below minimum wavelength")),
361  schema.addField<double>("throughputAtMax", "throughput above minimum wavelength")) {
362  Impl::setupPersistence(schema, arrays);
363  schema.getCitizen().markPersistent();
364  }
365  };
366 
367  void write(OutputArchiveHandle& handle) const override {
368  auto const& keys = PersistenceHelper::get();
369  auto catalog = handle.makeCatalog(keys.schema);
370  auto record = catalog.addNew();
371  record->set(keys.throughputAtMin, _atBounds.first);
372  record->set(keys.throughputAtMax, _atBounds.second);
373  _impl.persist(*record, keys.arrays);
374  handle.saveCatalog(catalog);
375  }
376 
377  class Factory : public table::io::PersistableFactory {
378  public:
379  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
380  CatalogVector const& catalogs) const override {
381  auto const& keys = PersistenceHelper::get();
382  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
383  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
384  auto& record = catalogs.front().front();
385  return std::make_shared<InterpolatedTransmissionCurve>(
386  Impl::unpersist(record, keys.arrays),
387  std::make_pair(record.get(keys.throughputAtMin), record.get(keys.throughputAtMax)));
388  }
389 
390  Factory() : table::io::PersistableFactory(Impl::NAME) {}
391  };
392 
393  static Factory registration;
394 
395 private:
396  std::pair<double, double> _atBounds;
397  Impl _impl;
398 };
399 
400 template <typename Impl>
401 typename InterpolatedTransmissionCurve<Impl>::Factory InterpolatedTransmissionCurve<Impl>::registration;
402 
403 template class InterpolatedTransmissionCurve<Impl1d>;
404 template class InterpolatedTransmissionCurve<Impl2d>;
405 
406 /*
407  * ProductTransmissionCurve: default for TransmissionCurve::multipliedBy().
408  *
409  * This is a straightforward lazy-evaluation object. Its only state is the
410  * two operands it delegates to.
411  */
412 class ProductTransmissionCurve : public TransmissionCurve {
413 public:
414  static constexpr char const* NAME = "ProductTransmissionCurve";
415 
416  ProductTransmissionCurve(std::shared_ptr<TransmissionCurve const> a,
418  : _a(std::move(a)), _b(std::move(b)) {}
419 
420  std::pair<double, double> getWavelengthBounds() const override {
421  auto aWavelengthBounds = _a->getWavelengthBounds();
422  auto bWavelengthBounds = _b->getWavelengthBounds();
423  auto aThroughputAtBounds = _a->getThroughputAtBounds();
424  auto bThroughputAtBounds = _b->getThroughputAtBounds();
425 
426  auto determineWavelengthBound = [](double aWavelength, double bWavelength, double aThroughput,
427  double bThroughput, auto isFirstOuter) -> double {
428  // Use the outermost wavelength bound only if its throughput
429  // values are not being multiplied by zeros from the operand with
430  // the innermost wavelength bound.
431  if (isFirstOuter(aWavelength, bWavelength)) {
432  return (bThroughput == 0.0) ? bWavelength : aWavelength;
433  } else {
434  return (aThroughput == 0.0) ? aWavelength : bWavelength;
435  }
436  };
437 
438  return std::make_pair(determineWavelengthBound(aWavelengthBounds.first, bWavelengthBounds.first,
439  aThroughputAtBounds.first, bThroughputAtBounds.first,
441  determineWavelengthBound(aWavelengthBounds.second, bWavelengthBounds.second,
442  aThroughputAtBounds.second, bThroughputAtBounds.second,
444  }
445 
446  std::pair<double, double> getThroughputAtBounds() const override {
447  auto aAtBounds = _a->getThroughputAtBounds();
448  auto bAtBounds = _b->getThroughputAtBounds();
449  return std::make_pair(aAtBounds.first * bAtBounds.first, aAtBounds.second * bAtBounds.second);
450  }
451 
452  void sampleAt(lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths,
453  ndarray::Array<double, 1, 1> const& out) const override {
454  _a->sampleAt(position, wavelengths, out);
455  ndarray::Array<double, 1, 1> tmp = ndarray::allocate(wavelengths.getSize<0>());
456  _b->sampleAt(position, wavelengths, tmp);
457  out.deep() *= tmp;
458  }
459 
460  bool isPersistable() const noexcept override { return _a->isPersistable() && _b->isPersistable(); }
461 
462 protected:
463  std::string getPersistenceName() const override { return NAME; }
464 
465  struct PersistenceHelper {
466  table::Schema schema;
467  table::Key<int> a;
468  table::Key<int> b;
469 
470  static PersistenceHelper const& get() {
471  static PersistenceHelper const instance;
472  return instance;
473  }
474 
475  private:
476  PersistenceHelper()
477  : schema(),
478  a(schema.addField<int>("a", "archive ID of first operand")),
479  b(schema.addField<int>("b", "archive ID of second operand")) {
480  schema.getCitizen().markPersistent();
481  }
482  };
483 
484  class Factory : public table::io::PersistableFactory {
485  public:
486  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
487  CatalogVector const& catalogs) const override {
488  auto const& keys = PersistenceHelper::get();
489  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
490  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
491  auto const& record = catalogs.front().front();
492  return std::make_shared<ProductTransmissionCurve>(
493  archive.get<TransmissionCurve>(record.get(keys.a)),
494  archive.get<TransmissionCurve>(record.get(keys.b)));
495  }
496 
497  Factory() : table::io::PersistableFactory(NAME) {}
498  };
499 
500  void write(OutputArchiveHandle& handle) const override {
501  auto const& keys = PersistenceHelper::get();
502  auto catalog = handle.makeCatalog(keys.schema);
503  auto record = catalog.addNew();
504  record->set(keys.a, handle.put(_a));
505  record->set(keys.b, handle.put(_b));
506  handle.saveCatalog(catalog);
507  }
508 
509  static Factory registration;
510 
511 private:
514 };
515 
516 ProductTransmissionCurve::Factory ProductTransmissionCurve::registration;
517 
518 /*
519  * TransformedTransmissionCurve: default for TransmissionCurve::transform.
520  *
521  * This is a another straightforward lazy-evaluation object. Its only state
522  * is the two operands it delegates to.
523  */
524 class TransformedTransmissionCurve : public TransmissionCurve {
525 public:
526  static constexpr char const* NAME = "TransformedTransmissionCurve";
527 
528  TransformedTransmissionCurve(std::shared_ptr<TransmissionCurve const> nested,
530  : _nested(std::move(nested)), _transform(std::move(transform)) {}
531 
532  std::pair<double, double> getWavelengthBounds() const override { return _nested->getWavelengthBounds(); }
533 
534  std::pair<double, double> getThroughputAtBounds() const override {
535  return _nested->getThroughputAtBounds();
536  }
537 
538  void sampleAt(lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths,
539  ndarray::Array<double, 1, 1> const& out) const override {
540  return _nested->sampleAt(_transform->applyInverse(position), wavelengths, out);
541  }
542 
543  bool isPersistable() const noexcept override {
544  return _nested->isPersistable() && _transform->isPersistable();
545  }
546 
547 protected:
548  // transforming a TransformedTransmissionCurve composes the transforms
549  std::shared_ptr<TransmissionCurve const> _transformedByImpl(
550  std::shared_ptr<geom::TransformPoint2ToPoint2> transform) const override {
551  return std::make_shared<TransformedTransmissionCurve>(_nested, transform->then(*_transform));
552  }
553 
554  std::string getPersistenceName() const override { return NAME; }
555 
556  struct PersistenceHelper {
557  table::Schema schema;
558  table::Key<int> nested;
559  table::Key<int> transform;
560 
561  static PersistenceHelper const& get() {
562  static PersistenceHelper const instance;
563  return instance;
564  }
565 
566  private:
567  PersistenceHelper()
568  : schema(),
569  nested(schema.addField<int>("nested", "archive ID of the nested TransmissionCurve")),
570  transform(schema.addField<int>("transform", "archive ID of the coordinate transform")) {
571  schema.getCitizen().markPersistent();
572  }
573  };
574 
575  void write(OutputArchiveHandle& handle) const override {
576  auto const& keys = PersistenceHelper::get();
577  auto catalog = handle.makeCatalog(keys.schema);
578  auto record = catalog.addNew();
579  record->set(keys.nested, handle.put(_nested));
580  record->set(keys.transform, handle.put(_transform));
581  handle.saveCatalog(catalog);
582  }
583 
584  class Factory : public table::io::PersistableFactory {
585  public:
586  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
587  CatalogVector const& catalogs) const override {
588  auto const& keys = PersistenceHelper::get();
589  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
590  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
591  auto const& record = catalogs.front().front();
592  return std::make_shared<TransformedTransmissionCurve>(
593  archive.get<TransmissionCurve>(record.get(keys.nested)),
594  archive.get<geom::TransformPoint2ToPoint2>(record.get(keys.transform)));
595  }
596 
597  Factory() : table::io::PersistableFactory(NAME) {}
598  };
599 
600  static Factory registration;
601 
602 private:
605 };
606 
607 TransformedTransmissionCurve::Factory TransformedTransmissionCurve::registration;
608 
609 } // namespace
610 
611 /*
612  * TransmissionCurve itself
613  */
614 
615 std::shared_ptr<TransmissionCurve const> TransmissionCurve::makeIdentity() {
616  return IdentityTransmissionCurve::get();
617 }
618 
619 std::shared_ptr<TransmissionCurve const> TransmissionCurve::makeSpatiallyConstant(
620  ndarray::Array<double const, 1> const& throughput, ndarray::Array<double const, 1> const& wavelengths,
621  double throughputAtMin, double throughputAtMax) {
622  ::gsl_set_error_handler_off();
623  LSST_THROW_IF_NE(wavelengths.getSize<0>(), throughput.getSize<0>(), pex::exceptions::LengthError,
624  "Length of wavelength array (%d) does not match size of throughput array (%d)");
625  return std::make_shared<InterpolatedTransmissionCurve<Impl1d>>(
626  Impl1d(ndarray::copy(throughput), ndarray::copy(wavelengths)),
627  std::make_pair(throughputAtMin, throughputAtMax));
628 }
629 
630 std::shared_ptr<TransmissionCurve const> TransmissionCurve::makeRadial(
631  ndarray::Array<double const, 2> const& throughput, ndarray::Array<double const, 1> const& wavelengths,
632  ndarray::Array<double const, 1> const& radii, double throughputAtMin, double throughputAtMax) {
633  ::gsl_set_error_handler_off();
635  wavelengths.getSize<0>(), throughput.getSize<0>(), pex::exceptions::LengthError,
636  "Length of wavelength array (%d) does not match first dimension of of throughput array (%d)");
638  radii.getSize<0>(), throughput.getSize<1>(), pex::exceptions::LengthError,
639  "Length of radii array (%d) does not match second dimension of of throughput array (%d)");
640  // GSL wants a column major array (Array<T,2,-2>). But ndarray can only flatten row-major arrays
641  // (Array<T,2,2>). So we allocate a row-major array, assign the caller's throughput array to a
642  // transposed view of it, and then flatten the row-major array.
643  ndarray::Array<double, 2, 2> throughputTransposed = ndarray::allocate(throughput.getShape().reverse());
644  throughputTransposed.transpose() = throughput;
645  ndarray::Array<double, 1, 1> throughputFlat = ndarray::flatten<1>(throughputTransposed);
646  return std::make_shared<InterpolatedTransmissionCurve<Impl2d>>(
647  Impl2d(throughputFlat, ndarray::copy(wavelengths), ndarray::copy(radii)),
648  std::make_pair(throughputAtMin, throughputAtMax));
649 }
650 
651 std::shared_ptr<TransmissionCurve const> TransmissionCurve::multipliedBy(
652  TransmissionCurve const& other) const {
653  auto a = shared_from_this();
654  auto b = other.shared_from_this();
655  auto result = a->_multipliedByImpl(b);
656  if (result == nullptr) {
657  result = b->_multipliedByImpl(a);
658  }
659  if (result == nullptr) {
660  result = std::make_shared<ProductTransmissionCurve>(std::move(a), std::move(b));
661  }
662  return result;
663 }
664 
665 std::shared_ptr<TransmissionCurve const> TransmissionCurve::transformedBy(
667  return _transformedByImpl(std::move(transform));
668 }
669 
670 ndarray::Array<double, 1, 1> TransmissionCurve::sampleAt(
671  lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths) const {
672  ndarray::Array<double, 1, 1> out = ndarray::allocate(wavelengths.getSize<0>());
673  sampleAt(position, wavelengths, out);
674  return out;
675 }
676 
677 std::shared_ptr<TransmissionCurve const> TransmissionCurve::_transformedByImpl(
679  return std::make_shared<TransformedTransmissionCurve>(shared_from_this(), std::move(transform));
680 }
681 
682 std::shared_ptr<TransmissionCurve const> TransmissionCurve::_multipliedByImpl(
684  return nullptr;
685 }
686 
687 std::string TransmissionCurve::getPythonModule() const { return "lsst.afw.image"; }
688 
689 } // namespace image
690 } // namespace afw
691 } // namespace lsst
#define LSST_CHECK_GSL(type, status)
A spatially-varying transmission curve as a function of wavelength.
#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
table::Key< int > b
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
py::object result
Definition: schema.cc:418
int y
Definition: SpanSet.cc:49
table::Schema schema
table::Key< int > a
T free(T... args)
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
Definition: Transform.h:300
STL class.
T min(T... args)
A base class for image defects.
T make_pair(T... args)
T isfinite(T... args)
T infinity(T... args)
T max(T... args)
T move(T... args)
table::Key< double > throughputAtMin
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
STL class.
STL class.
STL class.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
EigenVector const & asEigen() const noexcept(IS_ELEMENT_NOTHROW_COPYABLE)
Return a fixed-size Eigen representation of the coordinate object.
ItemVariant const * other
Definition: Schema.cc:56
T transform(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
ArrayKeyVector arrays
table::Key< double > throughputAtMax
table::Key< int > nested