LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
ConstrainedPhotometryModel.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * This file is part of jointcal.
4 *
5 * Developed for the LSST Data Management System.
6 * This product includes software developed by the LSST Project
7 * (https://www.lsst.org).
8 * See the COPYRIGHT file at the top-level directory of this distribution
9 * for details of code ownership.
10 *
11 * This program is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program. If not, see <https://www.gnu.org/licenses/>.
23 */
24
25#include <map>
26#include <limits>
27#include <vector>
28#include <string>
29
30#include "lsst/log/Log.h"
31
32#include "astshim.h"
33#include "astshim/ChebyMap.h"
34#include "lsst/geom.h"
41
42namespace lsst {
43namespace jointcal {
44
46 Eigen::Index firstIndex) {
47 Eigen::Index index = firstIndex;
48 if (whatToFit.find("Model") == std::string::npos) {
49 LOGLS_WARN(_log, "assignIndices was called and Model is *not* in whatToFit");
50 return index;
51 }
52
53 // If we got here, "Model" is definitely in whatToFit.
54 _fittingChips = (whatToFit.find("ModelChip") != std::string::npos);
55 _fittingVisits = (whatToFit.find("ModelVisit") != std::string::npos);
56 // If nothing more than "Model" is specified, it means fit everything.
57 if ((!_fittingChips) && (!_fittingVisits)) {
58 _fittingChips = _fittingVisits = true;
59 }
60
61 if (_fittingChips) {
62 for (auto &idMapping : _chipMap) {
63 auto mapping = idMapping.second.get();
64 // Don't assign indices for fixed parameters.
65 if (mapping->isFixed()) continue;
66 mapping->setIndex(index);
67 index += mapping->getNpar();
68 }
69 }
70 if (_fittingVisits) {
71 for (auto &idMapping : _visitMap) {
72 auto mapping = idMapping.second.get();
73 mapping->setIndex(index);
74 index += mapping->getNpar();
75 }
76 }
77 for (auto &idMapping : _chipVisitMap) {
78 idMapping.second->setWhatToFit(_fittingChips, _fittingVisits);
79 }
80 return index;
81}
82
83void ConstrainedPhotometryModel::offsetParams(Eigen::VectorXd const &delta) {
84 if (_fittingChips) {
85 for (auto &idMapping : _chipMap) {
86 auto mapping = idMapping.second.get();
87 // Don't offset indices for fixed parameters.
88 if (mapping->isFixed()) continue;
89 mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
90 }
91 }
92 if (_fittingVisits) {
93 for (auto &idMapping : _visitMap) {
94 auto mapping = idMapping.second.get();
95 mapping->offsetParams(delta.segment(mapping->getIndex(), mapping->getNpar()));
96 }
97 }
98}
99
101 for (auto &idMapping : _chipMap) {
102 idMapping.second->freezeErrorTransform();
103 }
104 for (auto &idMapping : _visitMap) {
105 idMapping.second->freezeErrorTransform();
106 }
107}
108
110 auto mapping = findMapping(ccdImage);
111 mapping->getMappingIndices(indices);
112}
113
115 std::size_t total = 0;
116 for (auto &idMapping : _chipMap) {
117 total += idMapping.second->getNpar();
118 }
119 for (auto &idMapping : _visitMap) {
120 total += idMapping.second->getNpar();
121 }
122 return total;
123}
124
126 CcdImage const &ccdImage,
127 Eigen::VectorXd &derivatives) const {
128 auto mapping = findMapping(ccdImage);
129 mapping->computeParameterDerivatives(measuredStar, measuredStar.getInstFlux(), derivatives);
130}
131
132namespace {
133// Convert photoTransform's way of storing Chebyshev coefficients into the format wanted by ChebyMap.
134ndarray::Array<double, 2, 2> toChebyMapCoeffs(const std::shared_ptr<PhotometryTransformChebyshev>& transform) {
135 auto coeffs = transform->getCoefficients();
136 // 4 x nPar: ChebyMap wants rows that look like (a_ij, 1, i, j) for out += a_ij*T_i(x)*T_j(y)
137 ndarray::Array<double, 2, 2> chebyCoeffs =
138 allocate(ndarray::makeVector(transform->getNpar(), std::size_t(4)));
139 Eigen::VectorXd::Index k = 0;
140 auto order = transform->getOrder();
141 for (ndarray::Size j = 0; j <= order; ++j) {
142 ndarray::Size const iMax = order - j; // to save re-computing `i+j <= order` every inner step.
143 for (ndarray::Size i = 0; i <= iMax; ++i, ++k) {
144 chebyCoeffs[k][0] = coeffs[j][i];
145 chebyCoeffs[k][1] = 1;
146 chebyCoeffs[k][2] = i;
147 chebyCoeffs[k][3] = j;
148 }
149 }
150 return chebyCoeffs;
151}
152} // namespace
153
155 for (auto &idMapping : _chipMap) {
156 out << "Sensor: " << idMapping.first << std::endl;
157 idMapping.second->print(out);
158 out << std::endl;
159 }
160 out << std::endl;
161 for (auto &idMapping : _visitMap) {
162 out << "Visit: " << idMapping.first << std::endl;
163 idMapping.second->print(out);
164 out << std::endl;
165 }
166}
167
169 auto idMapping = _chipVisitMap.find(ccdImage.getHashKey());
170 if (idMapping == _chipVisitMap.end())
172 "ConstrainedPhotometryModel cannot find CcdImage " + ccdImage.getName());
173 return idMapping->second.get();
174}
175
176template <class ChipTransform, class VisitTransform, class ChipVisitMapping>
178 geom::Box2D const &focalPlaneBBox, int visitOrder) {
179 // keep track of which chip we want to constrain (the one closest to the middle of the focal plane)
180 double minRadius2 = std::numeric_limits<double>::infinity();
181 CcdIdType constrainedChip = -1;
182
183 // First initialize all visit and ccd transforms, before we make the ccdImage mappings.
184 for (auto const &ccdImage : ccdImageList) {
185 auto visit = ccdImage->getVisit();
186 auto chip = ccdImage->getCcdId();
187 auto visitPair = _visitMap.find(visit);
188 auto chipPair = _chipMap.find(chip);
189
190 // If the chip is not in the map, add it, otherwise continue.
191 if (chipPair == _chipMap.end()) {
192 auto center = ccdImage->getDetector()->getCenter(afw::cameraGeom::FOCAL_PLANE);
193 double radius2 = std::pow(center.getX(), 2) + std::pow(center.getY(), 2);
194 if (radius2 < minRadius2) {
195 minRadius2 = radius2;
196 constrainedChip = chip;
197 }
198 auto photoCalib = ccdImage->getPhotoCalib();
199 // Use the single-frame processing calibration from the PhotoCalib as the default.
200 auto chipTransform = std::make_unique<ChipTransform>(initialChipCalibration(photoCalib));
201 _chipMap[chip] = std::make_shared<PhotometryMapping>(std::move(chipTransform));
202 }
203 // If the visit is not in the map, add it, otherwise continue.
204 if (visitPair == _visitMap.end()) {
205 auto visitTransform = std::make_unique<VisitTransform>(visitOrder, focalPlaneBBox);
206 _visitMap[visit] = std::make_shared<PhotometryMapping>(std::move(visitTransform));
207 }
208 }
209
210 // Fix one chip mapping, to remove the degeneracy from the system.
211 _chipMap.at(constrainedChip)->setFixed(true);
212
213 // Now create the ccdImage mappings, which are combinations of the chip/visit mappings above.
214 for (auto const &ccdImage : ccdImageList) {
215 auto visit = ccdImage->getVisit();
216 auto chip = ccdImage->getCcdId();
217 _chipVisitMap.emplace(ccdImage->getHashKey(),
218 std::make_unique<ChipVisitMapping>(_chipMap[chip], _visitMap[visit]));
219 }
220 LOGLS_INFO(_log, "Got " << _chipMap.size() << " chip mappings and " << _visitMap.size()
221 << " visit mappings; holding chip " << constrainedChip << " fixed ("
222 << getTotalParameters() << " total parameters).");
223 LOGLS_DEBUG(_log, "CcdImage map has " << _chipVisitMap.size() << " mappings, with "
224 << _chipVisitMap.bucket_count() << " buckets and a load factor of "
226}
227
229 CcdImage const &ccdImage) const {
230 auto detector = ccdImage.getDetector();
231 auto ccdBBox = detector->getBBox();
232 auto *mapping = dynamic_cast<ChipVisitPhotometryMapping *>(findMapping(ccdImage));
233
234 // There should be no way in which we can get to this point and not have a ChipVisitMapping,
235 // so blow up if we don't.
236 assert(mapping != nullptr);
237 // We know it's a Chebyshev transform because we created it as such, so blow up if it's not.
238 auto visitPhotometryTransform = std::dynamic_pointer_cast<PhotometryTransformChebyshev>(
239 mapping->getVisitMapping()->getTransform());
240 assert(visitPhotometryTransform != nullptr);
241 auto focalBBox = visitPhotometryTransform->getBBox();
242
243 // Unravel our chebyshev coefficients to build an astshim::ChebyMap.
244 auto coeff_f = toChebyMapCoeffs(std::dynamic_pointer_cast<PhotometryTransformChebyshev>(
245 mapping->getVisitMapping()->getTransform()));
246 // Bounds are the bbox
247 std::vector<double> lowerBound = {focalBBox.getMinX(), focalBBox.getMinY()};
248 std::vector<double> upperBound = {focalBBox.getMaxX(), focalBBox.getMaxY()};
249 afw::geom::TransformPoint2ToGeneric visitTransform(ast::ChebyMap(coeff_f, 1, lowerBound, upperBound));
250
251 double chipConstant = mapping->getChipMapping()->getParameters()[0];
252
253 // Compute a box that covers the area of the ccd in focal plane coordinates.
254 // This is the box over which we want to compute the mean of the visit transform.
255 auto pixToFocal = detector->getTransform(afw::cameraGeom::PIXELS, afw::cameraGeom::FOCAL_PLANE);
256 geom::Box2D ccdBBoxInFocal;
257 for (auto const &point : pixToFocal->applyForward(geom::Box2D(ccdBBox).getCorners())) {
258 ccdBBoxInFocal.include(point);
259 }
260 double visitMean = visitPhotometryTransform->mean(ccdBBoxInFocal);
261
262 return {chipConstant, visitTransform, pixToFocal, visitMean};
263}
264
265// ConstrainedFluxModel methods
266
268 MeasuredStar const &measuredStar) const {
269 return transform(ccdImage, measuredStar) - measuredStar.getFittedStar()->getFlux();
270}
271
272double ConstrainedFluxModel::transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const {
273 auto mapping = findMapping(ccdImage);
274 return mapping->transform(measuredStar, measuredStar.getInstFlux());
275}
276
278 MeasuredStar const &measuredStar) const {
279 auto mapping = findMapping(ccdImage);
280 double tempErr = tweakFluxError(measuredStar);
281 return mapping->transformError(measuredStar, measuredStar.getInstFlux(), tempErr);
282}
283
285 auto ccdBBox = ccdImage.getDetector()->getBBox();
286 auto prep = prepPhotoCalib(ccdImage);
287
288 // The chip part is easy: zoom map with the single value as the "zoom" factor
290 ast::ZoomMap(1, prep.chipConstant));
291
292 // Now stitch them all together.
293 auto transform = prep.pixToFocal->then(prep.visitTransform)->then(zoomTransform);
294
295 // NOTE: TransformBoundedField does not implement mean(), so we have to compute it here.
296 double mean = prep.chipConstant * prep.visitMean;
297
298 auto boundedField = std::make_shared<afw::math::TransformBoundedField>(ccdBBox, *transform);
299 return std::make_shared<afw::image::PhotoCalib>(mean, ccdImage.getPhotoCalib()->getCalibrationErr(),
300 boundedField, false);
301}
302
304 out << "ConstrainedFluxModel:" << std::endl;
306}
307
308// ConstrainedMagnitudeModel methods
309
311 MeasuredStar const &measuredStar) const {
312 return transform(ccdImage, measuredStar) - measuredStar.getFittedStar()->getMag();
313}
314
316 MeasuredStar const &measuredStar) const {
317 auto mapping = findMapping(ccdImage);
318 return mapping->transform(measuredStar, measuredStar.getInstMag());
319}
320
322 MeasuredStar const &measuredStar) const {
323 auto mapping = findMapping(ccdImage);
324 double tempErr = tweakFluxError(measuredStar);
325 return mapping->transformError(measuredStar, measuredStar.getInstFlux(), tempErr);
326}
327
329 CcdImage const &ccdImage) const {
330 auto ccdBBox = ccdImage.getDetector()->getBBox();
331 auto prep = prepPhotoCalib(ccdImage);
332
333 using namespace std::string_literals; // for operator""s to convert string literal->std::string
335 ast::MathMap(1, 1, {"y=pow(10.0,x/-2.5)"s}, {"x=-2.5*log10(y)"s}));
336
337 // The chip part is easy: zoom map with the value (converted to a flux) as the "zoom" factor.
338 double chipCalibration = utils::ABMagnitudeToNanojansky(prep.chipConstant);
340 ast::ZoomMap(1, chipCalibration));
341
342 // Now stitch them all together.
343 auto transform = prep.pixToFocal->then(prep.visitTransform)->then(logTransform)->then(zoomTransform);
344
345 // NOTE: TransformBoundedField does not implement mean(), so we have to compute it here.
346 double mean = chipCalibration * std::pow(10, prep.visitMean / -2.5);
347
348 auto boundedField = std::make_shared<afw::math::TransformBoundedField>(ccdBBox, *transform);
349 return std::make_shared<afw::image::PhotoCalib>(mean, ccdImage.getPhotoCalib()->getCalibrationErr(),
350 boundedField, false);
351}
352
354 out << "ConstrainedMagnitudeModel (" << _chipVisitMap.size() << " composite mappings; " << _chipMap.size()
355 << " sensor mappings, " << _visitMap.size() << " visit mappings):" << std::endl;
357}
358
359// explicit instantiation of templated function, so pybind11 can
362 geom::Box2D const &, int);
365 CcdImageList const &, geom::Box2D const &, int);
366
367} // namespace jointcal
368} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition Log.h:659
#define LOGLS_INFO(logger, message)
Log a info-level message using an iostream-based interface.
Definition Log.h:639
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
Definition Log.h:619
Implementation of the Photometric Calibration class.
table::Key< int > transform
T at(T... args)
T bucket_count(T... args)
A ChebyMap is a form of Mapping which performs a Chebyshev polynomial transformation.
Definition ChebyMap.h:97
A MathMap is a Mapping which allows you to specify a set of forward and/or inverse transformation fun...
Definition MathMap.h:61
A Mapping which "zooms" a set of points about the origin by multiplying all coordinate values by the ...
Definition ZoomMap.h:45
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint,...
Definition Transform.h:68
A floating-point coordinate rectangle geometry.
Definition Box.h:413
void include(Point2D const &point) noexcept
Expand this to ensure that this->contains(point).
Definition Box.cc:380
Handler of an actual image from a single CCD.
Definition CcdImage.h:64
std::shared_ptr< afw::cameraGeom::Detector > getDetector() const
Definition CcdImage.h:150
std::shared_ptr< afw::image::PhotoCalib > getPhotoCalib() const
Return the exposure's photometric calibration.
Definition CcdImage.h:160
std::string getName() const
Return the _name that identifies this ccdImage.
Definition CcdImage.h:79
CcdImageKey getHashKey() const
Definition CcdImage.h:152
A two-level photometric transform: one for the ccd and one for the visit.
std::shared_ptr< afw::image::PhotoCalib > toPhotoCalib(CcdImage const &ccdImage) const override
Return the mapping of ccdImage represented as a PhotoCalib.
double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override
Compute the residual between the model applied to a star and its associated fittedStar.
double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override
Return the on-sky transformed flux uncertainty for measuredStar on ccdImage.
double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override
Return the on-sky transformed flux for measuredStar on ccdImage.
void print(std::ostream &out) const override
Print a string representation of the contents of this mapping, for debugging.
double computeResidual(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override
Compute the residual between the model applied to a star and its associated fittedStar.
void print(std::ostream &out) const override
Print a string representation of the contents of this mapping, for debugging.
double transformError(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override
Return the on-sky transformed flux uncertainty for measuredStar on ccdImage.
double transform(CcdImage const &ccdImage, MeasuredStar const &measuredStar) const override
Return the on-sky transformed flux for measuredStar on ccdImage.
std::shared_ptr< afw::image::PhotoCalib > toPhotoCalib(CcdImage const &ccdImage) const override
Return the mapping of ccdImage represented as a PhotoCalib.
PrepPhotoCalib prepPhotoCalib(CcdImage const &ccdImage) const
Helper for preparing toPhotoCalib()
void initialize(CcdImageList const &ccdImageList, geom::Box2D const &focalPlaneBBox, int visitOrder)
Initialize the chip, visit, and chipVisit mappings by creating appropriate transforms and mappings.
PhotometryMappingBase * findMapping(CcdImage const &ccdImage) const override
Return a pointer to the mapping associated with this ccdImage.
void freezeErrorTransform() override
Once this routine has been called, the error transform is not modified by offsetParams().
void offsetParams(Eigen::VectorXd const &delta) override
Offset the parameters by the provided amounts (by -delta).
void print(std::ostream &out) const override
Print a string representation of the contents of this mapping, for debugging.
std::size_t getTotalParameters() const override
Return the total number of parameters in this model.
Eigen::Index assignIndices(std::string const &whatToFit, Eigen::Index firstIndex) override
Assign indices in the full matrix to the parameters being fit in the mappings, starting at firstIndex...
void getMappingIndices(CcdImage const &ccdImage, IndexVector &indices) const override
Get how this set of parameters (of length Npar()) map into the "grand" fit.
void computeParameterDerivatives(MeasuredStar const &measuredStar, CcdImage const &ccdImage, Eigen::VectorXd &derivatives) const override
Compute the parametric derivatives of this model.
virtual double initialChipCalibration(std::shared_ptr< afw::image::PhotoCalib const > photoCalib)=0
Return the initial calibration to use from this photoCalib.
nth-order 2d Chebyshev photometry transform, times the input flux.
Photometric offset independent of position, defined as (fluxMag0)^-1.
nth-order 2d Chebyshev photometry transform, plus the input flux.
Photometric offset independent of position, defined as -2.5 * log(flux / fluxMag0).
Sources measured on images.
std::shared_ptr< FittedStar > getFittedStar() const
Relates transform(s) to their position in the fitting matrix and allows interaction with the transfor...
LOG_LOGGER _log
lsst.logging instance, to be created by a subclass so that messages have consistent name.
double tweakFluxError(jointcal::MeasuredStar const &measuredStar) const
Add a fraction of the instrumental flux to the instrumental flux error, in quadrature.
Reports invalid arguments.
Definition Runtime.h:66
T emplace(T... args)
T end(T... args)
T endl(T... args)
T find(T... args)
T infinity(T... args)
T load_factor(T... args)
T move(T... args)
double ABMagnitudeToNanojansky(double magnitude)
Convert an AB magnitude to a flux in nanojansky.
Definition Magnitude.cc:32
T pow(T... args)
T size(T... args)
table::Key< int > order