LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
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 
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  }
364  };
365 
366  void write(OutputArchiveHandle& handle) const override {
367  auto const& keys = PersistenceHelper::get();
368  auto catalog = handle.makeCatalog(keys.schema);
369  auto record = catalog.addNew();
370  record->set(keys.throughputAtMin, _atBounds.first);
371  record->set(keys.throughputAtMax, _atBounds.second);
372  _impl.persist(*record, keys.arrays);
373  handle.saveCatalog(catalog);
374  }
375 
376  class Factory : public table::io::PersistableFactory {
377  public:
378  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
379  CatalogVector const& catalogs) const override {
380  auto const& keys = PersistenceHelper::get();
381  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
382  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
383  auto& record = catalogs.front().front();
384  return std::make_shared<InterpolatedTransmissionCurve>(
385  Impl::unpersist(record, keys.arrays),
386  std::make_pair(record.get(keys.throughputAtMin), record.get(keys.throughputAtMax)));
387  }
388 
389  Factory() : table::io::PersistableFactory(Impl::NAME) {}
390  };
391 
392  static Factory registration;
393 
394 private:
395  std::pair<double, double> _atBounds;
396  Impl _impl;
397 };
398 
399 template <typename Impl>
400 typename InterpolatedTransmissionCurve<Impl>::Factory InterpolatedTransmissionCurve<Impl>::registration;
401 
402 template class InterpolatedTransmissionCurve<Impl1d>;
403 template class InterpolatedTransmissionCurve<Impl2d>;
404 
405 /*
406  * ProductTransmissionCurve: default for TransmissionCurve::multipliedBy().
407  *
408  * This is a straightforward lazy-evaluation object. Its only state is the
409  * two operands it delegates to.
410  */
411 class ProductTransmissionCurve : public TransmissionCurve {
412 public:
413  static constexpr char const* NAME = "ProductTransmissionCurve";
414 
415  ProductTransmissionCurve(std::shared_ptr<TransmissionCurve const> a,
417  : _a(std::move(a)), _b(std::move(b)) {}
418 
419  std::pair<double, double> getWavelengthBounds() const override {
420  auto aWavelengthBounds = _a->getWavelengthBounds();
421  auto bWavelengthBounds = _b->getWavelengthBounds();
422  auto aThroughputAtBounds = _a->getThroughputAtBounds();
423  auto bThroughputAtBounds = _b->getThroughputAtBounds();
424 
425  auto determineWavelengthBound = [](double aWavelength, double bWavelength, double aThroughput,
426  double bThroughput, auto isFirstOuter) -> double {
427  // Use the outermost wavelength bound only if its throughput
428  // values are not being multiplied by zeros from the operand with
429  // the innermost wavelength bound.
430  if (isFirstOuter(aWavelength, bWavelength)) {
431  return (bThroughput == 0.0) ? bWavelength : aWavelength;
432  } else {
433  return (aThroughput == 0.0) ? aWavelength : bWavelength;
434  }
435  };
436 
437  return std::make_pair(determineWavelengthBound(aWavelengthBounds.first, bWavelengthBounds.first,
438  aThroughputAtBounds.first, bThroughputAtBounds.first,
440  determineWavelengthBound(aWavelengthBounds.second, bWavelengthBounds.second,
441  aThroughputAtBounds.second, bThroughputAtBounds.second,
443  }
444 
445  std::pair<double, double> getThroughputAtBounds() const override {
446  auto aAtBounds = _a->getThroughputAtBounds();
447  auto bAtBounds = _b->getThroughputAtBounds();
448  return std::make_pair(aAtBounds.first * bAtBounds.first, aAtBounds.second * bAtBounds.second);
449  }
450 
451  void sampleAt(lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths,
452  ndarray::Array<double, 1, 1> const& out) const override {
453  _a->sampleAt(position, wavelengths, out);
454  ndarray::Array<double, 1, 1> tmp = ndarray::allocate(wavelengths.getSize<0>());
455  _b->sampleAt(position, wavelengths, tmp);
456  out.deep() *= tmp;
457  }
458 
459  bool isPersistable() const noexcept override { return _a->isPersistable() && _b->isPersistable(); }
460 
461 protected:
462  std::string getPersistenceName() const override { return NAME; }
463 
464  struct PersistenceHelper {
465  table::Schema schema;
466  table::Key<int> a;
467  table::Key<int> b;
468 
469  static PersistenceHelper const& get() {
470  static PersistenceHelper const instance;
471  return instance;
472  }
473 
474  private:
475  PersistenceHelper()
476  : schema(),
477  a(schema.addField<int>("a", "archive ID of first operand")),
478  b(schema.addField<int>("b", "archive ID of second operand")) {}
479  };
480 
481  class Factory : public table::io::PersistableFactory {
482  public:
483  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
484  CatalogVector const& catalogs) const override {
485  auto const& keys = PersistenceHelper::get();
486  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
487  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
488  auto const& record = catalogs.front().front();
489  return std::make_shared<ProductTransmissionCurve>(
490  archive.get<TransmissionCurve>(record.get(keys.a)),
491  archive.get<TransmissionCurve>(record.get(keys.b)));
492  }
493 
494  Factory() : table::io::PersistableFactory(NAME) {}
495  };
496 
497  void write(OutputArchiveHandle& handle) const override {
498  auto const& keys = PersistenceHelper::get();
499  auto catalog = handle.makeCatalog(keys.schema);
500  auto record = catalog.addNew();
501  record->set(keys.a, handle.put(_a));
502  record->set(keys.b, handle.put(_b));
503  handle.saveCatalog(catalog);
504  }
505 
506  static Factory registration;
507 
508 private:
511 };
512 
513 ProductTransmissionCurve::Factory ProductTransmissionCurve::registration;
514 
515 /*
516  * TransformedTransmissionCurve: default for TransmissionCurve::transform.
517  *
518  * This is a another straightforward lazy-evaluation object. Its only state
519  * is the two operands it delegates to.
520  */
521 class TransformedTransmissionCurve : public TransmissionCurve {
522 public:
523  static constexpr char const* NAME = "TransformedTransmissionCurve";
524 
525  TransformedTransmissionCurve(std::shared_ptr<TransmissionCurve const> nested,
527  : _nested(std::move(nested)), _transform(std::move(transform)) {}
528 
529  std::pair<double, double> getWavelengthBounds() const override { return _nested->getWavelengthBounds(); }
530 
531  std::pair<double, double> getThroughputAtBounds() const override {
532  return _nested->getThroughputAtBounds();
533  }
534 
535  void sampleAt(lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths,
536  ndarray::Array<double, 1, 1> const& out) const override {
537  return _nested->sampleAt(_transform->applyInverse(position), wavelengths, out);
538  }
539 
540  bool isPersistable() const noexcept override {
541  return _nested->isPersistable() && _transform->isPersistable();
542  }
543 
544 protected:
545  // transforming a TransformedTransmissionCurve composes the transforms
546  std::shared_ptr<TransmissionCurve const> _transformedByImpl(
547  std::shared_ptr<geom::TransformPoint2ToPoint2> transform) const override {
548  return std::make_shared<TransformedTransmissionCurve>(_nested, transform->then(*_transform));
549  }
550 
551  std::string getPersistenceName() const override { return NAME; }
552 
553  struct PersistenceHelper {
554  table::Schema schema;
555  table::Key<int> nested;
556  table::Key<int> transform;
557 
558  static PersistenceHelper const& get() {
559  static PersistenceHelper const instance;
560  return instance;
561  }
562 
563  private:
564  PersistenceHelper()
565  : schema(),
566  nested(schema.addField<int>("nested", "archive ID of the nested TransmissionCurve")),
567  transform(schema.addField<int>("transform", "archive ID of the coordinate transform")) {}
568  };
569 
570  void write(OutputArchiveHandle& handle) const override {
571  auto const& keys = PersistenceHelper::get();
572  auto catalog = handle.makeCatalog(keys.schema);
573  auto record = catalog.addNew();
574  record->set(keys.nested, handle.put(_nested));
575  record->set(keys.transform, handle.put(_transform));
576  handle.saveCatalog(catalog);
577  }
578 
579  class Factory : public table::io::PersistableFactory {
580  public:
581  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
582  CatalogVector const& catalogs) const override {
583  auto const& keys = PersistenceHelper::get();
584  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
585  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
586  auto const& record = catalogs.front().front();
587  return std::make_shared<TransformedTransmissionCurve>(
588  archive.get<TransmissionCurve>(record.get(keys.nested)),
589  archive.get<geom::TransformPoint2ToPoint2>(record.get(keys.transform)));
590  }
591 
592  Factory() : table::io::PersistableFactory(NAME) {}
593  };
594 
595  static Factory registration;
596 
597 private:
600 };
601 
602 TransformedTransmissionCurve::Factory TransformedTransmissionCurve::registration;
603 
604 } // namespace
605 
606 /*
607  * TransmissionCurve itself
608  */
609 
610 std::shared_ptr<TransmissionCurve const> TransmissionCurve::makeIdentity() {
611  return IdentityTransmissionCurve::get();
612 }
613 
614 std::shared_ptr<TransmissionCurve const> TransmissionCurve::makeSpatiallyConstant(
615  ndarray::Array<double const, 1> const& throughput, ndarray::Array<double const, 1> const& wavelengths,
616  double throughputAtMin, double throughputAtMax) {
617  ::gsl_set_error_handler_off();
618  LSST_THROW_IF_NE(wavelengths.getSize<0>(), throughput.getSize<0>(), pex::exceptions::LengthError,
619  "Length of wavelength array (%d) does not match size of throughput array (%d)");
620  return std::make_shared<InterpolatedTransmissionCurve<Impl1d>>(
621  Impl1d(ndarray::copy(throughput), ndarray::copy(wavelengths)),
622  std::make_pair(throughputAtMin, throughputAtMax));
623 }
624 
625 std::shared_ptr<TransmissionCurve const> TransmissionCurve::makeRadial(
626  ndarray::Array<double const, 2> const& throughput, ndarray::Array<double const, 1> const& wavelengths,
627  ndarray::Array<double const, 1> const& radii, double throughputAtMin, double throughputAtMax) {
628  ::gsl_set_error_handler_off();
630  wavelengths.getSize<0>(), throughput.getSize<0>(), pex::exceptions::LengthError,
631  "Length of wavelength array (%d) does not match first dimension of of throughput array (%d)");
633  radii.getSize<0>(), throughput.getSize<1>(), pex::exceptions::LengthError,
634  "Length of radii array (%d) does not match second dimension of of throughput array (%d)");
635  // GSL wants a column major array (Array<T,2,-2>). But ndarray can only flatten row-major arrays
636  // (Array<T,2,2>). So we allocate a row-major array, assign the caller's throughput array to a
637  // transposed view of it, and then flatten the row-major array.
638  ndarray::Array<double, 2, 2> throughputTransposed = ndarray::allocate(throughput.getShape().reverse());
639  throughputTransposed.transpose() = throughput;
640  ndarray::Array<double, 1, 1> throughputFlat = ndarray::flatten<1>(throughputTransposed);
641  return std::make_shared<InterpolatedTransmissionCurve<Impl2d>>(
642  Impl2d(throughputFlat, ndarray::copy(wavelengths), ndarray::copy(radii)),
643  std::make_pair(throughputAtMin, throughputAtMax));
644 }
645 
646 std::shared_ptr<TransmissionCurve const> TransmissionCurve::multipliedBy(
647  TransmissionCurve const& other) const {
648  auto a = shared_from_this();
649  auto b = other.shared_from_this();
650  auto result = a->_multipliedByImpl(b);
651  if (result == nullptr) {
652  result = b->_multipliedByImpl(a);
653  }
654  if (result == nullptr) {
655  result = std::make_shared<ProductTransmissionCurve>(std::move(a), std::move(b));
656  }
657  return result;
658 }
659 
660 std::shared_ptr<TransmissionCurve const> TransmissionCurve::transformedBy(
662  return _transformedByImpl(std::move(transform));
663 }
664 
665 ndarray::Array<double, 1, 1> TransmissionCurve::sampleAt(
666  lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths) const {
667  ndarray::Array<double, 1, 1> out = ndarray::allocate(wavelengths.getSize<0>());
668  sampleAt(position, wavelengths, out);
669  return out;
670 }
671 
672 std::shared_ptr<TransmissionCurve const> TransmissionCurve::_transformedByImpl(
674  return std::make_shared<TransformedTransmissionCurve>(shared_from_this(), std::move(transform));
675 }
676 
677 std::shared_ptr<TransmissionCurve const> TransmissionCurve::_multipliedByImpl(
679  return nullptr;
680 }
681 
682 std::string TransmissionCurve::getPythonModule() const { return "lsst.afw.image"; }
683 
684 } // namespace image
685 } // namespace afw
686 } // namespace lsst
def write(self, patchRef, catalog)
Write the output.
#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
int y
Definition: SpanSet.cc:49
table::Schema schema
table::Key< int > a
ItemVariant const * other
Definition: Schema.cc:56
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.
T transform(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
ArrayKeyVector arrays
py::object result
Definition: _schema.cc:429
table::Key< double > throughputAtMax
table::Key< int > nested