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
PhotometryFit.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 <iostream>
26#include <iomanip>
27#include <algorithm>
28#include <fstream>
29#include <cmath>
30
31#include "Eigen/Sparse"
32
33#include "lsst/log/Log.h"
34#include "lsst/pex/exceptions.h"
37#include "lsst/jointcal/Chi2.h"
41
42namespace lsst {
43namespace jointcal {
44
46 Eigen::VectorXd &grad,
47 MeasuredStarList const *measuredStarList) const {
48 /**********************************************************************/
49 /* @note the math in this method and accumulateStatImageList() must be kept consistent,
50 * in terms of +/- convention, definition of model, etc. */
51 /**********************************************************************/
52
53 /* this routine works in two different ways: either providing the
54 Ccd, of providing the MeasuredStarList. In the latter case, the
55 Ccd should match the one(s) in the list. */
56 if (measuredStarList) assert(&(measuredStarList->front()->getCcdImage()) == &ccdImage);
57
58 std::size_t nparModel = (_fittingModel) ? _photometryModel->getNpar(ccdImage) : 0;
59 std::size_t nparFlux = (_fittingFluxes) ? 1 : 0;
60 std::size_t nparTotal = nparModel + nparFlux;
61 IndexVector indices(nparModel, -1);
62 if (_fittingModel) _photometryModel->getMappingIndices(ccdImage, indices);
63
64 Eigen::VectorXd H(nparTotal); // derivative matrix
65 // current position in the Jacobian
66 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
67 const MeasuredStarList &catalog = (measuredStarList) ? *measuredStarList : ccdImage.getCatalogForFit();
68
69 for (auto const &measuredStar : catalog) {
70 if (!measuredStar->isValid()) continue;
71 H.setZero(); // we cannot be sure that all entries will be overwritten.
72
73 double residual = _photometryModel->computeResidual(ccdImage, *measuredStar);
74 double inverseSigma = 1.0 / _photometryModel->transformError(ccdImage, *measuredStar);
75 double W = std::pow(inverseSigma, 2);
76
77 if (_fittingModel) {
78 _photometryModel->computeParameterDerivatives(*measuredStar, ccdImage, H);
79 for (std::size_t k = 0; k < indices.size(); k++) {
80 Eigen::Index l = indices[k];
81 tripletList.addTriplet(l, kTriplets, H[k] * inverseSigma);
82 grad[l] += H[k] * W * residual;
83 }
84 }
85 if (_fittingFluxes) {
86 Eigen::Index index = measuredStar->getFittedStar()->getIndexInMatrix();
87 // Note: H = dR/dFittedStarFlux == -1
88 tripletList.addTriplet(index, kTriplets, -1.0 * inverseSigma);
89 grad[index] += -1.0 * W * residual;
90 }
91 kTriplets += 1; // each measurement contributes 1 column in the Jacobian
92 }
93
94 tripletList.setNextFreeIndex(kTriplets);
95}
96
98 TripletList &tripletList, Eigen::VectorXd &grad) const {
99 /**********************************************************************/
102 /**********************************************************************/
103
104 // Derivatives of terms involving fitted and refstars only contribute if we are fitting fluxes.
105 if (!_fittingFluxes) return;
106 // Can't compute anything if there are no refStars.
107 if (_associations->refStarList.empty()) return;
108
109 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
110
111 for (auto const &fittedStar : fittedStarList) {
112 auto refStar = fittedStar->getRefStar();
113 if (refStar == nullptr) continue; // no contribution if no associated refstar
114 // TODO: Can we actually work with multiple filters at a time in this scheme?
115 // TODO: I feel we might only be able to fit one filter at a time, because these terms are
116 // independent of the ccdImage, so we don't have a specific filter.
117 // filter = ccdImage.getFilter();
118
119 // W == inverseSigma^2
120
121 double inverseSigma = 1.0 / _photometryModel->getRefError(*refStar);
122 // Residual is fittedStar - refStar for consistency with measurement terms.
123 double residual = _photometryModel->computeRefResidual(*fittedStar, *refStar);
124
125 Eigen::Index index = fittedStar->getIndexInMatrix();
126 // Note: H = dR/dFittedStar == 1
127 tripletList.addTriplet(index, kTriplets, 1.0 * inverseSigma);
128 grad(index) += 1.0 * std::pow(inverseSigma, 2) * residual;
129 kTriplets += 1;
130 }
131 tripletList.setNextFreeIndex(kTriplets);
132}
133
135 /**********************************************************************/
138 /**********************************************************************/
139 for (auto const &ccdImage : ccdImageList) {
140 auto &catalog = ccdImage->getCatalogForFit();
141
142 for (auto const &measuredStar : catalog) {
143 if (!measuredStar->isValid()) continue;
144 double sigma = _photometryModel->transformError(*ccdImage, *measuredStar);
145 double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
146
147 double chi2Val = std::pow(residual / sigma, 2);
148 accum.addEntry(chi2Val, 1, measuredStar);
149 } // end loop on measurements
150 }
151}
152
154 /**********************************************************************/
157 /**********************************************************************/
158
159 FittedStarList &fittedStarList = _associations->fittedStarList;
160 for (auto const &fittedStar : fittedStarList) {
161 auto refStar = fittedStar->getRefStar();
162 if (refStar == nullptr) continue;
163 double sigma = _photometryModel->getRefError(*refStar);
164 double residual = _photometryModel->computeRefResidual(*fittedStar, *refStar);
165 double chi2 = std::pow(residual / sigma, 2);
166 accum.addEntry(chi2, 1, fittedStar);
167 }
168}
169
171
173void PhotometryFit::getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const {
174 indices.clear();
175 if (_fittingModel) {
176 _photometryModel->getMappingIndices(measuredStar.getCcdImage(), indices);
177 }
178 if (_fittingFluxes) {
179 std::shared_ptr<FittedStar const> const fs = measuredStar.getFittedStar();
180 Eigen::Index fsIndex = fs->getIndexInMatrix();
181 indices.push_back(fsIndex);
182 }
183}
184
186 _whatToFit = whatToFit;
187 LOGLS_INFO(_log, "assignIndices: now fitting: " << whatToFit);
188 _fittingModel = (_whatToFit.find("Model") != std::string::npos);
189 _fittingFluxes = (_whatToFit.find("Fluxes") != std::string::npos);
190 // When entering here, we assume that whatToFit has already been interpreted.
191
192 _nModelParams = (_fittingModel) ? _photometryModel->assignIndices(whatToFit, 0) : 0;
194
195 if (_fittingFluxes) {
196 for (auto &fittedStar : _associations->fittedStarList) {
197 // the parameter layout here is used also
198 // - when filling the derivatives
199 // - when updating (offsetParams())
200 // - in getIndicesOfMeasuredStar
201 fittedStar->setIndexInMatrix(ipar);
202 ipar += 1;
203 }
204 }
206 _nTotal = ipar;
207 LOGLS_DEBUG(_log, "nParameters total: " << _nTotal << " model: " << _nModelParams
208 << " fluxes: " << _nStarParams);
209}
210
211void PhotometryFit::offsetParams(Eigen::VectorXd const &delta) {
212 if (delta.size() != _nTotal)
214 "PhotometryFit::offsetParams : the provided vector length is not compatible with "
215 "the current whatToFit setting");
216 if (_fittingModel) _photometryModel->offsetParams(delta);
217
218 if (_fittingFluxes) {
219 for (auto &fittedStar : _associations->fittedStarList) {
220 // the parameter layout here is used also
221 // - when filling the derivatives
222 // - when assigning indices (assignIndices())
223 Eigen::Index index = fittedStar->getIndexInMatrix();
224 _photometryModel->offsetFittedStar(*fittedStar, delta(index));
225 }
226 }
227}
228
230 std::ofstream ofile(filename.c_str());
231 std::string separator = "\t";
232
233 ofile << "id" << separator << "xccd" << separator << "yccd" << separator;
234 ofile << "mag" << separator << "instMag" << separator << "instMagErr" << separator;
235 ofile << "instFlux" << separator << "instFluxErr" << separator;
236 ofile << "inputFlux" << separator << "inputFluxErr" << separator;
237 ofile << "transformedFlux" << separator << "transformedFluxErr" << separator;
238 ofile << "fittedFlux" << separator;
239 ofile << "epoch" << separator;
240 ofile << "fsindex" << separator;
241 ofile << "ra" << separator << "dec" << separator;
242 ofile << "chi2" << separator << "nm" << separator;
243 ofile << "detector" << separator << "visit" << separator << std::endl;
244
245 ofile << "id in source catalog" << separator << "coordinates in CCD" << separator << separator;
246 ofile << "fitted magnitude" << separator << "measured magnitude" << separator
247 << "measured magnitude error" << separator;
248 ofile << "measured instrumental flux (ADU)" << separator << "measured instrument flux error" << separator;
249 ofile << "measured flux (nJy)" << separator << "measured flux error" << separator;
250 ofile << separator << separator;
251 ofile << "fitted flux (nJy)" << separator;
252 ofile << "Julian Epoch Year of the measurement" << separator;
253 ofile << "unique index of the fittedStar" << separator;
254 ofile << "on-sky position of fitted star" << separator << separator;
255 ofile << "contribution to chi2 (1 dof)" << separator << "number of measurements of this FittedStar"
256 << separator;
257 ofile << "detector id" << separator << "visit id" << std::endl;
258
259 const CcdImageList &ccdImageList = _associations->getCcdImageList();
260 for (auto const &ccdImage : ccdImageList) {
261 const MeasuredStarList &cat = ccdImage->getCatalogForFit();
262 for (auto const &measuredStar : cat) {
263 if (!measuredStar->isValid()) continue;
264
265 double instFluxErr = _photometryModel->tweakFluxError(*measuredStar);
266 double flux = _photometryModel->transform(*ccdImage, *measuredStar);
267 double fluxErr = _photometryModel->transformError(*ccdImage, *measuredStar);
268 std::shared_ptr<FittedStar const> const fittedStar = measuredStar->getFittedStar();
269 double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
270 double chi2Val = std::pow(residual / fluxErr, 2);
271
272 ofile << std::setprecision(9);
273 ofile << measuredStar->getId() << separator << measuredStar->x << separator << measuredStar->y
274 << separator;
275 ofile << fittedStar->getMag() << separator << measuredStar->getInstMag() << separator
276 << measuredStar->getInstMagErr() << separator;
277 ofile << measuredStar->getInstFlux() << separator << instFluxErr << separator;
278 ofile << measuredStar->getFlux() << separator << measuredStar->getFluxErr() << separator;
279 ofile << flux << separator << fluxErr << separator << fittedStar->getFlux() << separator;
280 ofile << ccdImage->getEpoch() << separator;
281 ofile << fittedStar->getIndexInMatrix() << separator;
282 ofile << fittedStar->x << separator << fittedStar->y << separator;
283 ofile << chi2Val << separator << fittedStar->getMeasurementCount() << separator;
284 ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() << std::endl;
285 } // loop on measurements in image
286 } // loop on images
287}
288
290 std::ofstream ofile(filename.c_str());
291 std::string separator = "\t";
292
293 ofile << "ra" << separator << "dec " << separator;
294 ofile << "mag" << separator;
295 ofile << "refFlux" << separator << "refFluxErr" << separator;
296 ofile << "fittedFlux" << separator << "fittedFluxErr" << separator;
297 ofile << "fsindex" << separator << "chi2" << separator << "nm" << std::endl;
298
299 ofile << "coordinates of fittedStar" << separator << separator;
300 ofile << "magnitude" << separator;
301 ofile << "refStar flux (nJy)" << separator << "refStar fluxErr" << separator;
302 ofile << "fittedStar flux (nJy)" << separator << "fittedStar fluxErr" << separator;
303 ofile << "unique index of the fittedStar" << separator << "refStar contribution to chi2 (1 dof)"
304 << separator << "number of measurements of this FittedStar" << std::endl;
305
306 // The following loop is heavily inspired from PhotometryFit::computeChi2()
307 const FittedStarList &fittedStarList = _associations->fittedStarList;
308 for (auto const &fittedStar : fittedStarList) {
309 const RefStar *refStar = fittedStar->getRefStar();
310 if (refStar == nullptr) continue;
311
312 double chi2 = std::pow(((fittedStar->getFlux() - refStar->getFlux()) / refStar->getFluxErr()), 2);
313
314 ofile << std::setprecision(9);
315 ofile << fittedStar->x << separator << fittedStar->y << separator;
316 ofile << fittedStar->getMag() << separator;
317 ofile << refStar->getFlux() << separator << refStar->getFluxErr() << separator;
318 ofile << fittedStar->getFlux() << separator << fittedStar->getFluxErr() << separator;
319 ofile << fittedStar->getIndexInMatrix() << separator << chi2 << separator
320 << fittedStar->getMeasurementCount() << std::endl;
321 } // loop on FittedStars
322}
323
324} // namespace jointcal
325} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< double > sigma
LSST DM logging module built on log4cxx.
#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
T c_str(T... args)
double getFlux() const
Definition BaseStar.h:98
double getFluxErr() const
Definition BaseStar.h:102
Handler of an actual image from a single CCD.
Definition CcdImage.h:64
MeasuredStarList const & getCatalogForFit() const
Gets the catalog to be used for fitting, which may have been cleaned-up.
Definition CcdImage.h:94
Base class for Chi2Statistic and Chi2List, to allow addEntry inside Fitter for either class.
Definition Chi2.h:44
virtual void addEntry(double inc, std::size_t dof, std::shared_ptr< BaseStar > star)=0
A list of FittedStar s. Such a list is typically constructed by Associations.
Definition FittedStar.h:116
std::shared_ptr< Associations > _associations
Definition FitterBase.h:165
Sources measured on images.
CcdImage const & getCcdImage() const
std::shared_ptr< FittedStar > getFittedStar() const
A list of MeasuredStar. They are usually filled in Associations::createCcdImage.
void accumulateStatRefStars(Chi2Accumulator &accum) const override
Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.
void assignIndices(std::string const &whatToFit) override
Set parameters to fit and assign indices in the big matrix.
void accumulateStatImageList(CcdImageList const &ccdImageList, Chi2Accumulator &accum) const override
Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for measurements.
void leastSquareDerivativesReference(FittedStarList const &fittedStarList, TripletList &tripletList, Eigen::VectorXd &grad) const override
Compute the derivatives of the reference terms.
void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const override
this routine is to be used only in the framework of outlier removal
void saveChi2RefContributions(std::string const &filename) const override
Save a CSV file containing residuals of reference terms.
void leastSquareDerivativesMeasurement(CcdImage const &ccdImage, TripletList &tripletList, Eigen::VectorXd &grad, MeasuredStarList const *measuredStarList=nullptr) const override
Compute the derivatives of the measured stars and model for one CcdImage.
void saveChi2MeasContributions(std::string const &filename) const override
Save a CSV file containing residuals of measurement terms.
void offsetParams(Eigen::VectorXd const &delta) override
Offset the parameters by the requested quantities.
Objects used as position/flux anchors (e.g.
Definition RefStar.h:45
Eigen::Index getNextFreeIndex() const
Definition Tripletlist.h:47
void setNextFreeIndex(Eigen::Index index)
Definition Tripletlist.h:49
void addTriplet(Eigen::Index i, Eigen::Index j, double val)
Definition Tripletlist.h:43
Reports invalid arguments.
Definition Runtime.h:66
T clear(T... args)
T endl(T... args)
T find(T... args)
T front(T... args)
T pow(T... args)
T push_back(T... args)
T setprecision(T... args)
T size(T... args)