LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
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 
42 namespace lsst {
43 namespace jointcal {
44 
45 void PhotometryFit::leastSquareDerivativesMeasurement(CcdImage const &ccdImage, TripletList &tripletList,
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 
97 void PhotometryFit::leastSquareDerivativesReference(FittedStarList const &fittedStarList,
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.size() == 0) 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 
134 void PhotometryFit::accumulateStatImageList(CcdImageList const &ccdImageList, Chi2Accumulator &accum) const {
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 
153 void PhotometryFit::accumulateStatRefStars(Chi2Accumulator &accum) const {
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 
173 void PhotometryFit::getIndicesOfMeasuredStar(MeasuredStar const &measuredStar,
174  IndexVector &indices) const {
175  indices.clear();
176  if (_fittingModel) {
177  _photometryModel->getMappingIndices(measuredStar.getCcdImage(), indices);
178  }
179  if (_fittingFluxes) {
180  std::shared_ptr<FittedStar const> const fs = measuredStar.getFittedStar();
181  Eigen::Index fsIndex = fs->getIndexInMatrix();
182  indices.push_back(fsIndex);
183  }
184 }
185 
187  _whatToFit = whatToFit;
188  LOGLS_INFO(_log, "assignIndices: now fitting: " << whatToFit);
189  _fittingModel = (_whatToFit.find("Model") != std::string::npos);
190  _fittingFluxes = (_whatToFit.find("Fluxes") != std::string::npos);
191  // When entering here, we assume that whatToFit has already been interpreted.
192 
193  _nParModel = (_fittingModel) ? _photometryModel->assignIndices(whatToFit, 0) : 0;
194  std::size_t ipar = _nParModel;
195 
196  if (_fittingFluxes) {
197  for (auto &fittedStar : _associations->fittedStarList) {
198  // the parameter layout here is used also
199  // - when filling the derivatives
200  // - when updating (offsetParams())
201  // - in getIndicesOfMeasuredStar
202  fittedStar->setIndexInMatrix(ipar);
203  ipar += 1;
204  }
205  }
206  _nParFluxes = ipar - _nParModel;
207  _nParTot = ipar;
209  "nParameters total: " << _nParTot << " model: " << _nParModel << " fluxes: " << _nParFluxes);
210 }
211 
212 void PhotometryFit::offsetParams(Eigen::VectorXd const &delta) {
213  if (delta.size() != _nParTot)
215  "PhotometryFit::offsetParams : the provided vector length is not compatible with "
216  "the current whatToFit setting");
217  if (_fittingModel) _photometryModel->offsetParams(delta);
218 
219  if (_fittingFluxes) {
220  for (auto &fittedStar : _associations->fittedStarList) {
221  // the parameter layout here is used also
222  // - when filling the derivatives
223  // - when assigning indices (assignIndices())
224  Eigen::Index index = fittedStar->getIndexInMatrix();
225  _photometryModel->offsetFittedStar(*fittedStar, delta(index));
226  }
227  }
228 }
229 
231  std::ofstream ofile(filename.c_str());
232  std::string separator = "\t";
233 
234  ofile << "#id" << separator << "xccd" << separator << "yccd" << separator;
235  ofile << "mag" << separator << "instMag" << separator << "instMagErr" << separator;
236  ofile << "instFlux" << separator << "instFluxErr" << separator;
237  ofile << "inputFlux" << separator << "inputFluxErr" << separator;
238  ofile << "transformedFlux" << separator << "transformedFluxErr" << separator;
239  ofile << "fittedFlux" << separator;
240  ofile << "mjd" << separator << "color" << separator;
241  ofile << "fsindex" << separator;
242  ofile << "ra" << separator << "dec" << separator;
243  ofile << "chi2" << separator << "nm" << separator;
244  ofile << "chip" << separator << "visit" << separator << std::endl;
245 
246  ofile << "#id in source catalog" << separator << "coordinates in CCD" << separator << separator;
247  ofile << "fitted magnitude" << separator << "measured magnitude" << separator
248  << "measured magnitude error" << separator;
249  ofile << "measured instrumental flux (ADU)" << separator << "measured instrument flux error" << separator;
250  ofile << "measured flux (nJy)" << separator << "measured flux error" << separator;
251  ofile << separator << separator;
252  ofile << "fitted flux (nJy)" << separator;
253  ofile << "modified Julian date of the measurement" << separator << "currently unused" << separator;
254  ofile << "unique index of the fittedStar" << separator;
255  ofile << "on-sky position of fitted star" << separator << separator;
256  ofile << "contribution to Chi2 (1 dof)" << separator << "number of measurements of this FittedStar"
257  << separator;
258  ofile << "chip id" << separator << "visit id" << std::endl;
259 
260  const CcdImageList &ccdImageList = _associations->getCcdImageList();
261  for (auto const &ccdImage : ccdImageList) {
262  const MeasuredStarList &cat = ccdImage->getCatalogForFit();
263  for (auto const &measuredStar : cat) {
264  if (!measuredStar->isValid()) continue;
265 
266  double instFluxErr = _photometryModel->tweakFluxError(*measuredStar);
267  double flux = _photometryModel->transform(*ccdImage, *measuredStar);
268  double fluxErr = _photometryModel->transformError(*ccdImage, *measuredStar);
269  double jd = ccdImage->getMjd();
270  std::shared_ptr<FittedStar const> const fittedStar = measuredStar->getFittedStar();
271  double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
272  double chi2Val = std::pow(residual / fluxErr, 2);
273 
274  ofile << std::setprecision(9);
275  ofile << measuredStar->getId() << separator << measuredStar->x << separator << measuredStar->y
276  << separator;
277  ofile << fittedStar->getMag() << separator << measuredStar->getInstMag() << separator
278  << measuredStar->getInstMagErr() << separator;
279  ofile << measuredStar->getInstFlux() << separator << instFluxErr << separator;
280  ofile << measuredStar->getFlux() << separator << measuredStar->getFluxErr() << separator;
281  ofile << flux << separator << fluxErr << separator << fittedStar->getFlux() << separator;
282  ofile << jd << separator << fittedStar->color << separator;
283  ofile << fittedStar->getIndexInMatrix() << separator;
284  ofile << fittedStar->x << separator << fittedStar->y << separator;
285  ofile << chi2Val << separator << fittedStar->getMeasurementCount() << separator;
286  ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() << std::endl;
287  } // loop on measurements in image
288  } // loop on images
289 }
290 
292  std::ofstream ofile(filename.c_str());
293  std::string separator = "\t";
294 
295  ofile << "#ra" << separator << "dec " << separator;
296  ofile << "mag" << separator << "color" << separator;
297  ofile << "refFlux" << separator << "refFluxErr" << separator;
298  ofile << "fittedFlux" << separator << "fittedFluxErr" << separator;
299  ofile << "fsindex" << separator << "chi2" << separator << "nm" << std::endl;
300 
301  ofile << "#coordinates of fittedStar" << separator << separator;
302  ofile << "magnitude" << separator << "currently unused" << separator;
303  ofile << "refStar flux (nJy)" << separator << "refStar fluxErr" << separator;
304  ofile << "fittedStar flux (nJy)" << separator << "fittedStar fluxErr" << separator;
305  ofile << "unique index of the fittedStar" << separator << "refStar contribution to Chi2 (1 dof)"
306  << separator << "number of measurements of this FittedStar" << std::endl;
307 
308  // The following loop is heavily inspired from PhotometryFit::computeChi2()
309  const FittedStarList &fittedStarList = _associations->fittedStarList;
310  for (auto const &fittedStar : fittedStarList) {
311  const RefStar *refStar = fittedStar->getRefStar();
312  if (refStar == nullptr) continue;
313 
314  double chi2 = std::pow(((fittedStar->getFlux() - refStar->getFlux()) / refStar->getFluxErr()), 2);
315 
316  ofile << std::setprecision(9);
317  ofile << fittedStar->x << separator << fittedStar->y << separator;
318  ofile << fittedStar->getMag() << separator << fittedStar->color << separator;
319  ofile << refStar->getFlux() << separator << refStar->getFluxErr() << separator;
320  ofile << fittedStar->getFlux() << separator << fittedStar->getFluxErr() << separator;
321  ofile << fittedStar->getIndexInMatrix() << separator << chi2 << separator
322  << fittedStar->getMeasurementCount() << std::endl;
323  } // loop on FittedStars
324 }
325 
326 } // namespace jointcal
327 } // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
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:628
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
Definition: Log.h:608
T c_str(T... args)
double getFlux() const
Definition: BaseStar.h:96
double getFluxErr() const
Definition: BaseStar.h:100
A list of FittedStar s. Such a list is typically constructed by Associations.
Definition: FittedStar.h:123
std::shared_ptr< Associations > _associations
Definition: FitterBase.h:153
A list of MeasuredStar. They are usually filled in Associations::createCcdImage.
Definition: MeasuredStar.h:146
void assignIndices(std::string const &whatToFit) override
Set parameters to fit and assign indices in the big matrix.
void saveChi2RefContributions(std::string const &filename) const override
Save a CSV file containing residuals of reference terms.
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 anchors, typically USNO stars.
Definition: RefStar.h:39
Reports invalid arguments.
Definition: Runtime.h:66
T clear(T... args)
T endl(T... args)
T find(T... args)
std::list< std::shared_ptr< CcdImage > > CcdImageList
Definition: CcdImage.h:46
A base class for image defects.
T pow(T... args)
T push_back(T... args)
T setprecision(T... args)