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
AstrometryFit.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 
30 #include "Eigen/Sparse"
31 
32 #include "lsst/log/Log.h"
33 #include "lsst/pex/exceptions.h"
36 #include "lsst/jointcal/Chi2.h"
42 
43 namespace {
53 double computeProjectedRefStarChi2(lsst::jointcal::FatPoint refStar) {
54  double det = refStar.vx * refStar.vy - std::pow(refStar.vxy, 2);
55  double wxx = refStar.vy / det;
56  double wyy = refStar.vx / det;
57  double wxy = -refStar.vxy / det;
58  return wxx * std::pow(refStar.x, 2) + 2 * wxy * refStar.x * refStar.y + wyy * std::pow(refStar.y, 2);
59 }
60 } // namespace
61 
62 namespace lsst {
63 namespace jointcal {
64 
66  std::shared_ptr<AstrometryModel> astrometryModel, double posError)
67  : FitterBase(associations),
68  _astrometryModel(astrometryModel),
69  _epoch(_associations->getEpoch()),
70  _posError(posError) {
71  _log = LOG_GET("lsst.jointcal.AstrometryFit");
72 }
73 
74 /* ! this routine is used in 3 instances: when computing
75 the derivatives, when computing the Chi2, when filling a tuple.
76 */
77 Point AstrometryFit::transformFittedStar(FittedStar const &fittedStar, AstrometryTransform const &sky2TP,
78  double deltaYears) const {
79  Point fittedStarInTP;
80  if (fittedStar.getRefStar()) {
81  // PM data is on-sky, not in the tangent plane
82  Point temp = fittedStar.getRefStar()->applyProperMotion(fittedStar, deltaYears);
83  fittedStarInTP = sky2TP.apply(temp);
84  } else {
85  fittedStarInTP = sky2TP.apply(fittedStar);
86  }
87  return fittedStarInTP;
88 }
89 
93 static void tweakAstromMeasurementErrors(FatPoint &P, MeasuredStar const &Ms, double error) {
94  static bool called = false;
95  static double increment = 0;
96  if (!called) {
97  increment = std::pow(error, 2); // was in Preferences
98  called = true;
99  }
100  P.vx += increment;
101  P.vy += increment;
102 }
103 
104 // we could consider computing the chi2 here.
105 // (although it is not extremely useful)
106 void AstrometryFit::leastSquareDerivativesMeasurement(CcdImage const &ccdImage, TripletList &tripletList,
107  Eigen::VectorXd &fullGrad,
108  MeasuredStarList const *msList) const {
109  /**********************************************************************/
110  /* @note the math in this method and accumulateStatImage() must be kept consistent,
111  * in terms of +/- convention, definition of model, etc. */
112  /**********************************************************************/
113 
114  /* Setup */
115  /* this routine works in two different ways: either providing the
116  ccdImage, of providing the MeasuredStarList. In the latter case, the
117  ccdImage should match the one(s) in the list. */
118  if (msList) assert(&(msList->front()->getCcdImage()) == &ccdImage);
119 
120  // get the Mapping
121  const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
122  // count parameters
123  std::size_t npar_mapping = (_fittingDistortions) ? mapping->getNpar() : 0;
124  std::size_t npar_pos = (_fittingPos) ? 2 : 0;
125  std::size_t npar_pm = (_fittingPM) ? 2 : 0;
126  std::size_t npar_tot = npar_mapping + npar_pos + npar_pm;
127  // if (npar_tot == 0) this CcdImage does not contribute
128  // any constraint to the fit, so :
129  if (npar_tot == 0) return;
130  IndexVector indices(npar_tot, -1);
131  if (_fittingDistortions) mapping->getMappingIndices(indices);
132 
133  // FittedStar is "observed" epoch, MeasuredStar is "baseline"
134  double deltaYears = _epoch - ccdImage.getEpoch();
135  // transformation from sky to TP
136  auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
137  // reserve matrices once for all measurements
138  AstrometryTransformLinear dypdy;
139  // the shape of H (et al) is required this way in order to be able to
140  // separate derivatives along x and y as vectors.
141  Eigen::MatrixX2d H(npar_tot, 2), halpha(npar_tot, 2), HW(npar_tot, 2);
142  Eigen::Matrix2d transW(2, 2);
143  Eigen::Matrix2d alpha(2, 2);
144  Eigen::VectorXd grad(npar_tot);
145  // current position in the Jacobian
146  Eigen::Index kTriplets = tripletList.getNextFreeIndex();
147  const MeasuredStarList &catalog = (msList) ? *msList : ccdImage.getCatalogForFit();
148 
149  for (auto &i : catalog) {
150  const MeasuredStar &ms = *i;
151  if (!ms.isValid()) continue;
152  // tweak the measurement errors
153  FatPoint inPos = ms;
154  tweakAstromMeasurementErrors(inPos, ms, _posError);
155  H.setZero(); // we cannot be sure that all entries will be overwritten.
156  FatPoint outPos;
157  // should *not* fill H if whatToFit excludes mapping parameters.
158  if (_fittingDistortions)
159  mapping->computeTransformAndDerivatives(inPos, outPos, H);
160  else
161  mapping->transformPosAndErrors(inPos, outPos);
162 
163  std::size_t ipar = npar_mapping;
164  double det = outPos.vx * outPos.vy - std::pow(outPos.vxy, 2);
165  if (det <= 0 || outPos.vx <= 0 || outPos.vy <= 0) {
166  LOGLS_WARN(_log, "Inconsistent measurement errors: drop measurement at "
167  << Point(ms) << " in image " << ccdImage.getName());
168  continue;
169  }
170  transW(0, 0) = outPos.vy / det;
171  transW(1, 1) = outPos.vx / det;
172  transW(0, 1) = transW(1, 0) = -outPos.vxy / det;
173  // compute alpha, a triangular square root
174  // of transW (i.e. a Cholesky factor)
175  alpha(0, 0) = sqrt(transW(0, 0));
176  // checked that alpha*alphaT = transW
177  alpha(1, 0) = transW(0, 1) / alpha(0, 0);
178  // DB - I think that the next line is equivalent to : alpha(1,1) = 1./sqrt(outPos.vy)
179  // PA - seems correct !
180  alpha(1, 1) = 1. / sqrt(det * transW(0, 0));
181  alpha(0, 1) = 0;
182 
183  std::shared_ptr<FittedStar const> const fs = ms.getFittedStar();
184 
185  Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
186 
187  // compute derivative of TP position w.r.t sky position ....
188  if (npar_pos > 0) // ... if actually fitting FittedStar position
189  {
190  sky2TP->computeDerivative(*fs, dypdy, 1e-3);
191  // sign checked
192  // TODO Still have to check with non trivial non-diagonal terms
193  H(npar_mapping, 0) = -dypdy.A11();
194  H(npar_mapping + 1, 0) = -dypdy.A12();
195  H(npar_mapping, 1) = -dypdy.A21();
196  H(npar_mapping + 1, 1) = -dypdy.A22();
197  indices[npar_mapping] = fs->getIndexInMatrix();
198  indices.at(npar_mapping + 1) = fs->getIndexInMatrix() + 1;
199  ipar += npar_pos;
200  }
201  /* only consider proper motions of objects allowed to move,
202  unless the fit is going to be degenerate */
203  // TODO: left as reference for when we implement PM fitting
204  // TODO: mjd here would become either deltaYears or maybe associations.epoch? Check the math first!
205  // if (_fittingPM && fs->mightMove) {
206  // H(ipar, 0) = -mjd; // Sign unchecked but consistent with above
207  // H(ipar + 1, 1) = -mjd;
208  // indices[ipar] = fs->getIndexInMatrix() + 2;
209  // indices[ipar + 1] = fs->getIndexInMatrix() + 3;
210  // ipar += npar_pm;
211  // }
212 
213  // We can now compute the residual
214  Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
215 
216  // do not write grad = H*transW*res to avoid
217  // dynamic allocation of a temporary
218  halpha = H * alpha;
219  HW = H * transW;
220  grad = HW * res;
221  // now feed in triplets and fullGrad
222  for (std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
223  for (std::size_t ic = 0; ic < 2; ++ic) {
224  double val = halpha(ipar, ic);
225  if (val == 0) continue;
226  tripletList.addTriplet(indices[ipar], kTriplets + ic, val);
227  }
228  fullGrad(indices[ipar]) += grad(ipar);
229  }
230  kTriplets += 2; // each measurement contributes 2 columns in the Jacobian
231  } // end loop on measurements
232  tripletList.setNextFreeIndex(kTriplets);
233 }
234 
235 void AstrometryFit::leastSquareDerivativesReference(FittedStarList const &fittedStarList,
236  TripletList &tripletList,
237  Eigen::VectorXd &fullGrad) const {
238  /**********************************************************************/
239  /* @note the math in this method and accumulateStatRefStars() must be kept consistent,
240  * in terms of +/- convention, definition of model, etc. */
241  /**********************************************************************/
242 
243  /* We compute here the derivatives of the terms involving fitted
244  stars and reference stars. They only provide contributions if we
245  are fitting positions: */
246  if (!_fittingPos) return;
247  /* the other case where the accumulation of derivatives stops
248  here is when there are no RefStars */
249  if (_associations->refStarList.size() == 0) return;
250  Eigen::Matrix2d W(2, 2);
251  Eigen::Matrix2d alpha(2, 2);
252  Eigen::Matrix2d H(2, 2), halpha(2, 2), HW(2, 2);
253  AstrometryTransformLinear der;
254  Eigen::Vector2d res, grad;
255  Eigen::Index indices[2];
256  Eigen::Index kTriplets = tripletList.getNextFreeIndex();
257  /* We cannot use the spherical coordinates directly to evaluate
258  Euclidean distances, we have to use a projector on some plane in
259  order to express least squares. Not projecting could lead to a
260  disaster around the poles or across alpha=0. So we need a
261  projector. We construct a projector and will change its
262  projection point at every object */
263  TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
264  for (auto const &i : fittedStarList) {
265  const FittedStar &fs = *i;
266  auto rs = fs.getRefStar();
267  if (rs == nullptr) continue;
268  proj.setTangentPoint(fs);
269  // fs projects to (0,0), no need to compute its transform.
270  FatPoint rsProj;
271  proj.transformPosAndErrors(*rs, rsProj);
272  // Compute the derivative of the projector to incorporate its effects on the errors.
273  proj.computeDerivative(fs, der, 1e-4);
274  // sign checked. TODO check that the off-diagonal terms are OK.
275  H(0, 0) = -der.A11();
276  H(1, 0) = -der.A12();
277  H(0, 1) = -der.A21();
278  H(1, 1) = -der.A22();
279  // TO DO : account for proper motions.
280  double det = rsProj.vx * rsProj.vy - std::pow(rsProj.vxy, 2);
281  if (rsProj.vx <= 0 || rsProj.vy <= 0 || det <= 0) {
282  LOGLS_WARN(_log, "RefStar error matrix not positive definite for: " << *rs);
283  continue;
284  }
285  W(0, 0) = rsProj.vy / det;
286  W(0, 1) = W(1, 0) = -rsProj.vxy / det;
287  W(1, 1) = rsProj.vx / det;
288  // compute alpha, a triangular square root
289  // of W (i.e. a Cholesky factor)
290  alpha(0, 0) = sqrt(W(0, 0));
291  // checked that alpha*alphaT = transW
292  alpha(1, 0) = W(0, 1) / alpha(0, 0);
293  alpha(1, 1) = 1. / sqrt(det * W(0, 0));
294  alpha(0, 1) = 0;
295  indices[0] = fs.getIndexInMatrix();
296  indices[1] = fs.getIndexInMatrix() + 1;
297  unsigned npar_tot = 2;
298  /* TODO: account here for proper motions in the reference
299  catalog. We can code the effect and set the value to 0. Most
300  (all?) catalogs do not even come with a reference epoch. Gaia
301  will change that. When refraction enters into the game, one should
302  pay attention to the orientation of the frame */
303 
304  /* The residual should be Proj(fs)-Proj(*rs) in order to be consistent
305  with the measurement terms. Since P(fs) = 0, we have: */
306  res[0] = -rsProj.x;
307  res[1] = -rsProj.y;
308  halpha = H * alpha;
309  // grad = H*W*res
310  HW = H * W;
311  grad = HW * res;
312  // now feed in triplets and fullGrad
313  for (std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
314  for (unsigned ic = 0; ic < 2; ++ic) {
315  double val = halpha(ipar, ic);
316  if (val == 0) continue;
317  tripletList.addTriplet(indices[ipar], kTriplets + ic, val);
318  }
319  fullGrad(indices[ipar]) += grad(ipar);
320  }
321  kTriplets += 2; // each measurement contributes 2 columns in the Jacobian
322  }
323  tripletList.setNextFreeIndex(kTriplets);
324 }
325 
326 void AstrometryFit::accumulateStatImage(CcdImage const &ccdImage, Chi2Accumulator &accum) const {
327  /**********************************************************************/
330  /**********************************************************************/
331  /* Setup */
332  // 1 : get the Mapping's
333  const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
334  // FittedStar is "observed" epoch, MeasuredStar is "baseline"
335  double deltaYears = _epoch - ccdImage.getEpoch();
336  // transformation from sky to TP
337  auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
338  // reserve matrix once for all measurements
339  Eigen::Matrix2Xd transW(2, 2);
340 
341  auto &catalog = ccdImage.getCatalogForFit();
342  for (auto const &ms : catalog) {
343  if (!ms->isValid()) continue;
344  // tweak the measurement errors
345  FatPoint inPos = *ms;
346  tweakAstromMeasurementErrors(inPos, *ms, _posError);
347 
348  FatPoint outPos;
349  // should *not* fill H if whatToFit excludes mapping parameters.
350  mapping->transformPosAndErrors(inPos, outPos);
351  double det = outPos.vx * outPos.vy - std::pow(outPos.vxy, 2);
352  if (det <= 0 || outPos.vx <= 0 || outPos.vy <= 0) {
353  LOGLS_WARN(_log, " Inconsistent measurement errors :drop measurement at "
354  << Point(*ms) << " in image " << ccdImage.getName());
355  continue;
356  }
357  transW(0, 0) = outPos.vy / det;
358  transW(1, 1) = outPos.vx / det;
359  transW(0, 1) = transW(1, 0) = -outPos.vxy / det;
360 
361  std::shared_ptr<FittedStar const> const fs = ms->getFittedStar();
362  Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
363 
364  Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
365  double chi2Val = res.transpose() * transW * res;
366 
367  accum.addEntry(chi2Val, 2, ms);
368  } // end of loop on measurements
369 }
370 
371 void AstrometryFit::accumulateStatImageList(CcdImageList const &ccdImageList, Chi2Accumulator &accum) const {
372  for (auto const &ccdImage : ccdImageList) {
373  accumulateStatImage(*ccdImage, accum);
374  }
375 }
376 
377 void AstrometryFit::accumulateStatRefStars(Chi2Accumulator &accum) const {
378  /**********************************************************************/
381  /**********************************************************************/
382 
383  /* If you wonder why we project here, read comments in
384  AstrometryFit::leastSquareDerivativesReference(TripletList &TList, Eigen::VectorXd &Rhs) */
385  FittedStarList &fittedStarList = _associations->fittedStarList;
386  TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
387  for (auto const &fs : fittedStarList) {
388  auto rs = fs->getRefStar();
389  if (rs == nullptr) continue;
390  proj.setTangentPoint(*fs);
391  // fs projects to (0,0), no need to compute its transform.
392  FatPoint rsProj;
393  proj.transformPosAndErrors(*rs, rsProj);
394  // TO DO : account for proper motions.
395  double chi2 = computeProjectedRefStarChi2(rsProj);
396  accum.addEntry(chi2, 2, fs);
397  }
398 }
399 
401 
403 void AstrometryFit::getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const {
404  if (_fittingDistortions) {
405  const AstrometryMapping *mapping = _astrometryModel->getMapping(measuredStar.getCcdImage());
406  mapping->getMappingIndices(indices);
407  }
408  std::shared_ptr<FittedStar const> const fs = measuredStar.getFittedStar();
409  Eigen::Index fsIndex = fs->getIndexInMatrix();
410  if (_fittingPos) {
411  indices.push_back(fsIndex);
412  indices.push_back(fsIndex + 1);
413  }
414  // For securing the outlier removal, the next block is just useless
415  if (_fittingPM) {
416  for (std::size_t k = 0; k < 2; ++k) indices.push_back(fsIndex + 2 + k);
417  }
418  /* Should not put the index of refaction stuff or we will not be
419  able to remove more than 1 star at a time. */
420 }
421 
423  _whatToFit = whatToFit;
424  LOGLS_INFO(_log, "assignIndices: Now fitting " << whatToFit);
425  _fittingDistortions = (_whatToFit.find("Distortions") != std::string::npos);
426  _fittingPos = (_whatToFit.find("Positions") != std::string::npos);
427  _fittingPM = (_whatToFit.find("PM") != std::string::npos);
428  // When entering here, we assume that whatToFit has already been interpreted.
429 
430  _nModelParams = 0;
431  if (_fittingDistortions) _nModelParams = _astrometryModel->assignIndices(_whatToFit, 0);
432  std::size_t ipar = _nModelParams;
433 
434  if (_fittingPos) {
435  FittedStarList &fittedStarList = _associations->fittedStarList;
436  for (auto &fittedStar : fittedStarList) {
437  // the parameter layout here is used also
438  // - when filling the derivatives
439  // - when updating (offsetParams())
440  // - in GetMeasuredStarIndices
441  fittedStar->setIndexInMatrix(ipar);
442  ipar += 2;
443  // TODO: left as reference for when we implement PM fitting
444  // if ((_fittingPM)&fittedStar->mightMove) ipar += 2;
445  }
446  }
447  _nStarParams = ipar - _nModelParams;
448  _nTotal = ipar;
449  LOGLS_DEBUG(_log, "nParameters total: " << _nTotal << " model: " << _nModelParams
450  << " values: " << _nStarParams);
451 }
452 
453 void AstrometryFit::offsetParams(Eigen::VectorXd const &delta) {
454  if (delta.size() != _nTotal)
456  "AstrometryFit::offsetParams : the provided vector length is not compatible with "
457  "the current whatToFit setting");
458  if (_fittingDistortions) _astrometryModel->offsetParams(delta);
459 
460  if (_fittingPos) {
461  FittedStarList &fittedStarList = _associations->fittedStarList;
462  for (auto const &i : fittedStarList) {
463  FittedStar &fs = *i;
464  // the parameter layout here is used also
465  // - when filling the derivatives
466  // - when assigning indices (assignIndices())
467  Eigen::Index index = fs.getIndexInMatrix();
468  fs.x += delta(index);
469  fs.y += delta(index + 1);
470  // TODO: left as reference for when we implement PM fitting
471  // if ((_fittingPM)&fs.mightMove) {
472  // fs.pmx += delta(index + 2);
473  // fs.pmy += delta(index + 3);
474  // }
475  }
476  }
477 }
478 
480  const char *what2fit[] = {"Positions", "Distortions", "Positions Distortions"};
481  // DEBUG
482  for (unsigned k = 0; k < sizeof(what2fit) / sizeof(what2fit[0]); ++k) {
483  assignIndices(what2fit[k]);
484  TripletList tripletList(10000);
485  Eigen::VectorXd grad(_nTotal);
486  grad.setZero();
487  leastSquareDerivatives(tripletList, grad);
488  SparseMatrixD jacobian(_nTotal, tripletList.getNextFreeIndex());
489  jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
490  SparseMatrixD hessian = jacobian * jacobian.transpose();
491  LOGLS_DEBUG(_log, "npar : " << _nTotal << ' ' << _nModelParams);
492  }
493 }
494 
496  std::ofstream ofile(filename.c_str());
497  std::string separator = "\t";
498 
499  ofile << "id" << separator << "xccd" << separator << "yccd " << separator;
500  ofile << "rx" << separator << "ry" << separator;
501  ofile << "xtp" << separator << "ytp" << separator;
502  ofile << "mag" << separator << "deltaYears" << separator;
503  ofile << "xErr" << separator << "yErr" << separator << "xyCov" << separator;
504  ofile << "xtpi" << separator << "ytpi" << separator;
505  ofile << "rxi" << separator << "ryi" << separator;
506  ofile << "fsindex" << separator;
507  ofile << "ra" << separator << "dec" << separator;
508  ofile << "chi2" << separator << "nm" << separator;
509  ofile << "detector" << separator << "visit" << std::endl;
510 
511  ofile << "id in source catalog" << separator << "coordinates in CCD (pixels)" << separator << separator;
512  ofile << "residual on TP (degrees)" << separator << separator;
513  ofile << "transformed coordinate in TP (degrees)" << separator << separator;
514  ofile << "rough magnitude" << separator << "Julian epoch year delta from fit epoch" << separator;
515  ofile << "transformed measurement uncertainty (degrees)" << separator << separator << separator;
516  ofile << "as-read position on TP (degrees)" << separator << separator;
517  ofile << "as-read residual on TP (degrees)" << separator << separator;
518  ofile << "unique index of the fittedStar" << separator;
519  ofile << "on sky position of fittedStar" << separator << separator;
520  ofile << "contribution to chi2 (2D dofs)" << separator << "number of measurements of this fittedStar"
521  << separator;
522  ofile << "detector id" << separator << "visit id" << std::endl;
523 
524  const CcdImageList &ccdImageList = _associations->getCcdImageList();
525  for (auto const &ccdImage : ccdImageList) {
526  const MeasuredStarList &cat = ccdImage->getCatalogForFit();
527  const AstrometryMapping *mapping = _astrometryModel->getMapping(*ccdImage);
528  const auto readTransform = ccdImage->getReadWcs();
529  // FittedStar is "observed" epoch, MeasuredStar is "baseline"
530  double deltaYears = _epoch - ccdImage->getEpoch();
531  for (auto const &ms : cat) {
532  if (!ms->isValid()) continue;
533  FatPoint tpPos;
534  FatPoint inPos = *ms;
535  tweakAstromMeasurementErrors(inPos, *ms, _posError);
536  mapping->transformPosAndErrors(inPos, tpPos);
537  auto sky2TP = _astrometryModel->getSkyToTangentPlane(*ccdImage);
538  const std::unique_ptr<AstrometryTransform> readPixToTangentPlane =
539  compose(*sky2TP, *readTransform);
540  FatPoint inputTpPos = readPixToTangentPlane->apply(inPos);
541  std::shared_ptr<FittedStar const> const fs = ms->getFittedStar();
542 
543  Point fittedStarInTP = transformFittedStar(*fs, *sky2TP, deltaYears);
544  Point res = tpPos - fittedStarInTP;
545  Point inputRes = inputTpPos - fittedStarInTP;
546  double det = tpPos.vx * tpPos.vy - std::pow(tpPos.vxy, 2);
547  double wxx = tpPos.vy / det;
548  double wyy = tpPos.vx / det;
549  double wxy = -tpPos.vxy / det;
550  double chi2 = wxx * res.x * res.x + wyy * res.y * res.y + 2 * wxy * res.x * res.y;
551 
552  ofile << std::setprecision(9);
553  ofile << ms->getId() << separator << ms->x << separator << ms->y << separator;
554  ofile << res.x << separator << res.y << separator;
555  ofile << tpPos.x << separator << tpPos.y << separator;
556  ofile << fs->getMag() << separator << deltaYears << separator;
557  ofile << tpPos.vx << separator << tpPos.vy << separator << tpPos.vxy << separator;
558  ofile << inputTpPos.x << separator << inputTpPos.y << separator;
559  ofile << inputRes.x << separator << inputRes.y << separator;
560  ofile << fs->getIndexInMatrix() << separator;
561  ofile << fs->x << separator << fs->y << separator;
562  ofile << chi2 << separator << fs->getMeasurementCount() << separator;
563  ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() << std::endl;
564  } // loop on measurements in image
565  } // loop on images
566 }
567 
569  std::ofstream ofile(filename.c_str());
570  std::string separator = "\t";
571 
572  ofile << "ra" << separator << "dec " << separator;
573  ofile << "rx" << separator << "ry" << separator;
574  ofile << "mag" << separator;
575  ofile << "xErr" << separator << "yErr" << separator << "xyCov" << separator;
576  ofile << "fsindex" << separator;
577  ofile << "chi2" << separator << "nm" << std::endl;
578 
579  ofile << "coordinates of fittedStar (degrees)" << separator << separator;
580  ofile << "residual on TP (degrees)" << separator << separator;
581  ofile << "magnitude" << separator;
582  ofile << "refStar transformed measurement uncertainty (degrees)" << separator << separator << separator;
583  ofile << "unique index of the fittedStar" << separator;
584  ofile << "refStar contribution to chi2 (2D dofs)" << separator
585  << "number of measurements of this FittedStar" << std::endl;
586 
587  // The following loop is heavily inspired from AstrometryFit::computeChi2()
588  const FittedStarList &fittedStarList = _associations->fittedStarList;
590  for (auto const &i : fittedStarList) {
591  const FittedStar &fs = *i;
592  auto rs = fs.getRefStar();
593  if (rs == nullptr) continue;
594  proj.setTangentPoint(fs);
595  // fs projects to (0,0), no need to compute its transform.
596  FatPoint rsProj;
597  proj.transformPosAndErrors(*rs, rsProj);
598  double chi2 = computeProjectedRefStarChi2(rsProj);
599 
600  ofile << std::setprecision(9);
601  ofile << fs.x << separator << fs.y << separator;
602  ofile << rsProj.x << separator << rsProj.y << separator;
603  ofile << fs.getMag() << separator;
604  ofile << rsProj.vx << separator << rsProj.vy << separator << rsProj.vxy << separator;
605  ofile << fs.getIndexInMatrix() << separator;
606  ofile << chi2 << separator << fs.getMeasurementCount() << std::endl;
607  } // loop on FittedStars
608 }
609 } // namespace jointcal
610 } // namespace lsst
Eigen::Matrix< double, Eigen::Dynamic, 2 > MatrixX2d
Definition: Eigenstuff.h:33
Eigen::SparseMatrix< double, 0, Eigen::Index > SparseMatrixD
Definition: Eigenstuff.h:35
#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 LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
#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 begin(T... args)
T c_str(T... args)
void checkStuff()
DEBUGGING routine.
void saveChi2RefContributions(std::string const &filename) const override
Save a CSV file containing residuals of reference terms.
void offsetParams(Eigen::VectorXd const &delta) override
Offset the parameters by the requested quantities.
AstrometryFit(std::shared_ptr< Associations > associations, std::shared_ptr< AstrometryModel > astrometryModel, double posError)
this is the only constructor
void assignIndices(std::string const &whatToFit) override
Set parameters to fit and assign indices in the big matrix.
void saveChi2MeasContributions(std::string const &filename) const override
Save a CSV file containing residuals of measurement terms.
virtual class needed in the abstraction of the distortion model
virtual void transformPosAndErrors(FatPoint const &where, FatPoint &outPoint) const =0
The same as above but without the parameter derivatives (used to evaluate chi^2)
a virtual (interface) class for geometric transformations.
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const =0
virtual void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, const double step=0.01) const
Computes the local Derivative of a transform, w.r.t.
implements the linear transformations (6 real coefficients).
double getMag() const
Definition: BaseStar.h:105
A Point with uncertainties.
Definition: FatPoint.h:34
FittedStars are objects whose position or flux is going to be fitted, and which come from the associa...
Definition: FittedStar.h:50
const RefStar * getRefStar() const
Get the astrometric reference star associated with this star.
Definition: FittedStar.h:105
void setIndexInMatrix(Eigen::Index const index)
index is a value that a fit can set and reread....
Definition: FittedStar.h:96
int getMeasurementCount() const
The number of MeasuredStars currently associated with this FittedStar.
Definition: FittedStar.h:88
Eigen::Index getIndexInMatrix() const
Definition: FittedStar.h:99
A list of FittedStar s. Such a list is typically constructed by Associations.
Definition: FittedStar.h:116
Base class for fitters.
Definition: FitterBase.h:53
void leastSquareDerivatives(TripletList &tripletList, Eigen::VectorXd &grad) const
Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting.
Definition: FitterBase.cc:344
Eigen::Index _nModelParams
Definition: FitterBase.h:168
std::shared_ptr< Associations > _associations
Definition: FitterBase.h:163
A list of MeasuredStar. They are usually filled in Associations::createCcdImage.
Definition: MeasuredStar.h:151
A point in a plane.
Definition: Point.h:37
double x
coordinate
Definition: Point.h:42
Point applyProperMotion(Point star, double timeDeltaYears) const
Apply proper motion correction to the input star, returning a star with PM-corrected coordinates and ...
Definition: RefStar.cc:34
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane)
void transformPosAndErrors(const FatPoint &in, FatPoint &out) const
transform with analytical derivatives
void setTangentPoint(Point const &tangentPoint)
Resets the projection (or tangent) point.
Eigen::Index getNextFreeIndex() const
Definition: Tripletlist.h:47
Reports invalid arguments.
Definition: Runtime.h:66
T end(T... args)
T endl(T... args)
T find(T... args)
std::list< std::shared_ptr< CcdImage > > CcdImageList
Definition: CcdImage.h:46
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
A base class for image defects.
T pow(T... args)
T push_back(T... args)
T setprecision(T... args)
T sqrt(T... args)
ImageT val
Definition: CR.cc:146