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