LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
ProductBoundedField.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 
24 #include <algorithm>
25 
26 #include "ndarray/eigen.h"
27 
28 #include "lsst/pex/exceptions.h"
34 
35 namespace lsst {
36 namespace afw {
37 
38 template std::shared_ptr<math::ProductBoundedField> table::io::PersistableFacade<
39  math::ProductBoundedField>::dynamicCast(std::shared_ptr<table::io::Persistable> const&);
40 
41 namespace math {
42 
43 namespace {
44 
45 lsst::geom::Box2I checkAndExtractBBox(std::vector<std::shared_ptr<BoundedField const>> const & factors) {
46  if (factors.size() < 1u) {
47  throw LSST_EXCEPT(
48  pex::exceptions::LengthError,
49  "ProductBoundedField requires at least one BoundedField factor."
50  );
51  }
52  auto iter = factors.begin();
53  auto bbox = (**iter).getBBox();
54  ++iter;
55  for (; iter != factors.end(); ++iter) {
56  if ((**iter).getBBox() != bbox) {
57  throw LSST_EXCEPT(
58  pex::exceptions::InvalidParameterError,
59  (boost::format("Inconsistency in ProductBoundedField bboxes: %s != %s") %
60  bbox % (**iter).getBBox()).str()
61  );
62  }
63  }
64  return bbox;
65 }
66 
67 } // anonymous
68 
70  BoundedField(checkAndExtractBBox(factors)), _factors(factors)
71 {}
72 
76 
77 double ProductBoundedField::evaluate(lsst::geom::Point2D const& position) const {
78  double product = 1.0;
79  for (auto const & field : _factors) {
80  product *= field->evaluate(position);
81  }
82  return product;
83 }
84 
85 ndarray::Array<double, 1, 1> ProductBoundedField::evaluate(
86  ndarray::Array<double const, 1> const& x,
87  ndarray::Array<double const, 1> const& y
88 ) const {
89  if (x.getShape() != y.getShape()) {
90  throw LSST_EXCEPT(
92  (boost::format("Inconsistent shapes: %s != %s") % x.getShape() % y.getShape()).str()
93  );
94  }
95  ndarray::Array<double, 1, 1> z = ndarray::allocate(x.getShape());
96  std::fill(z.begin(), z.end(), 1.0);
97  for (auto const & field : _factors) {
98  ndarray::asEigenArray(z) *= ndarray::asEigenArray(field->evaluate(x, y));
99  }
100  return z;
101 }
102 
103 // ------------------ persistence ---------------------------------------------------------------------------
104 
105 namespace {
106 
107 struct PersistenceHelper {
108  table::Schema schema;
109  table::Key<int> id;
110 
111  static PersistenceHelper const & get() {
112  static PersistenceHelper const instance;
113  return instance;
114  }
115 
116 private:
117 
118  PersistenceHelper() :
119  schema(),
120  id(schema.addField<int>("id", "Archive ID of a BoundedField factor."))
121  {}
122 
123  PersistenceHelper(PersistenceHelper const &) = delete;
124  PersistenceHelper(PersistenceHelper &&) = delete;
125  PersistenceHelper & operator=(PersistenceHelper const &) = delete;
126  PersistenceHelper & operator=(PersistenceHelper &&) = delete;
127 
128  ~PersistenceHelper() noexcept = default;
129 
130 };
131 
132 class ProductBoundedFieldFactory : public table::io::PersistableFactory {
133 public:
134  explicit ProductBoundedFieldFactory(std::string const& name)
135  : afw::table::io::PersistableFactory(name) {}
136 
137  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
138  CatalogVector const& catalogs) const override {
139  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
140  auto const & keys = PersistenceHelper::get();
141  auto const & cat = catalogs.front();
143  factors.reserve(cat.size());
144  for (auto const & record : cat) {
145  factors.push_back(archive.get<BoundedField>(record.get(keys.id)));
146  }
147  return std::make_shared<ProductBoundedField>(factors);
148  }
149 };
150 
151 std::string getProductBoundedFieldPersistenceName() { return "ProductBoundedField"; }
152 
153 ProductBoundedFieldFactory registration(getProductBoundedFieldPersistenceName());
154 
155 } // namespace
156 
157 bool ProductBoundedField::isPersistable() const noexcept {
158  return std::all_of(_factors.begin(), _factors.end(),
159  [](auto const & field) { return field->isPersistable(); });
160 }
161 
163  return getProductBoundedFieldPersistenceName();
164 }
165 
166 std::string ProductBoundedField::getPythonModule() const { return "lsst.afw.math"; }
167 
169  auto const & keys = PersistenceHelper::get();
170  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
171  catalog.reserve(_factors.size());
172  for (auto const & field : _factors) {
173  catalog.addNew()->set(keys.id, handle.put(field));
174  }
175  handle.saveCatalog(catalog);
176 }
177 
178 // ------------------ operators -----------------------------------------------------------------------------
179 
182  bool multiplied = false;
183  for (auto & field : factors) {
184  try {
185  field = (*field) * scale;
186  multiplied = true;
187  break;
188  } catch (pex::exceptions::LogicError &) {}
189  }
190  if (!multiplied) {
191  ndarray::Array<double, 2, 2> coefficients = ndarray::allocate(1, 1);
192  coefficients[0][0] = scale;
193  factors.push_back(std::make_shared<ChebyshevBoundedField>(getBBox(), coefficients));
194  }
195  return std::make_shared<ProductBoundedField>(factors);
196 }
197 
199  auto rhsCasted = dynamic_cast<ProductBoundedField const*>(&rhs);
200  if (!rhsCasted) return false;
201 
202  return (getBBox() == rhsCasted->getBBox()) &&
203  std::equal(_factors.begin(), _factors.end(),
204  rhsCasted->_factors.begin(), rhsCasted->_factors.end(),
205  [](auto const & a, auto const & b) { return *a == *b; });
206 }
207 
208 std::string ProductBoundedField::toString() const {
210  os << "ProductBoundedField([";
211  for (auto const & field : _factors) {
212  os << (*field) << ", ";
213  }
214  os << "])";
215  return os.str();
216 }
217 
218 } // namespace math
219 } // namespace afw
220 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
table::Key< int > field
Definition: ApCorrMap.cc:77
ndarray::Array< double const, 2, 2 > coefficients
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
double z
Definition: Match.cc:44
table::Key< int > id
table::Schema schema
std::ostream * os
Definition: Schema.cc:557
int y
Definition: SpanSet.cc:48
table::Key< int > b
table::Key< int > a
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
T all_of(T... args)
An abstract base class for 2-d functions defined on an integer bounding boxes.
Definition: BoundedField.h:55
lsst::geom::Box2I getBBox() const
Return the bounding box that defines the region where the field is valid.
Definition: BoundedField.h:112
A BoundedField that lazily multiplies a sequence of other BoundedFields.
bool operator==(BoundedField const &rhs) const override
BoundedFields (of the same sublcass) are equal if their bounding boxes and parameters are equal.
bool isPersistable() const noexcept override
ProductBoundedField is persistable if and only if all of its factors are.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
ProductBoundedField(std::vector< std::shared_ptr< BoundedField const >> const &factors)
Construct from a sequence of BoundedField factors.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
std::shared_ptr< BoundedField > operator*(double const scale) const override
Return a scaled BoundedField.
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
virtual double evaluate(lsst::geom::Point2D const &position) const=0
Evaluate the field at the given point.
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:490
void reserve(size_type n)
Increase the capacity of the catalog to the given size.
Definition: Catalog.h:433
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
An integer coordinate rectangle.
Definition: Box.h:55
Reports invalid arguments.
Definition: Runtime.h:66
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
T equal(T... args)
T fill(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
FilterProperty & operator=(FilterProperty const &)=default
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
Definition: warpExposure.cc:0