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
FitterBase.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 <vector>
26 #include "Eigen/Core"
27 
28 #include <boost/math/tools/minima.hpp>
29 
30 #include "lsst/log/Log.h"
31 
32 #include "lsst/jointcal/Chi2.h"
33 #include "lsst/jointcal/CcdImage.h"
38 
39 namespace lsst {
40 namespace jointcal {
41 
43  Chi2Statistic chi2;
44  accumulateStatImageList(_associations->getCcdImageList(), chi2);
46  // chi2.ndof contains the number of squares.
47  // So subtract the number of parameters.
48  chi2.ndof -= _nTotal;
49  return chi2;
50 }
51 
53  FittedStarList &fsOutliers, double &cut) const {
54  // collect chi2 contributions
55  Chi2List chi2List;
56  chi2List.reserve(_associations->getMaxMeasuredStars() + _associations->refStarList.size());
57  // contributions from measurement terms:
58  accumulateStatImageList(_associations->ccdImageList, chi2List);
59  // and from reference terms
60  accumulateStatRefStars(chi2List);
61 
62  // compute some statistics
63  size_t nval = chi2List.size();
64  if (nval == 0) return 0;
65  sort(chi2List.begin(), chi2List.end());
66  double median = (nval & 1) ? chi2List[nval / 2].chi2
67  : 0.5 * (chi2List[nval / 2 - 1].chi2 + chi2List[nval / 2].chi2);
68  auto averageAndSigma = chi2List.computeAverageAndSigma();
69  LOGLS_DEBUG(_log, "findOutliers chi2 stat: mean/median/sigma " << averageAndSigma.first << '/' << median
70  << '/' << averageAndSigma.second);
71  cut = averageAndSigma.first + nSigmaCut * averageAndSigma.second;
72  /* For each of the parameters, we will not remove more than 1
73  measurement that contributes to constraining it. Keep track using
74  of what we are touching using an integer vector. This is the
75  trick that Marc Betoule came up to for outlier removals in "star
76  flats" fits. */
77  Eigen::VectorXi affectedParams(_nTotal);
78  affectedParams.setZero();
79 
80  std::size_t nOutliers = 0; // returned to the caller
81  // start from the strongest outliers.
82  for (auto chi2 = chi2List.rbegin(); chi2 != chi2List.rend(); ++chi2) {
83  if (chi2->chi2 < cut) break; // because the array is sorted.
84  IndexVector indices;
85  /* now, we want to get the indices of the parameters this chi2
86  term depends on. We have to figure out which kind of term it
87  is; we use for that the type of the star attached to the Chi2Star. */
88  auto measuredStar = std::dynamic_pointer_cast<MeasuredStar>(chi2->star);
89  std::shared_ptr<FittedStar> fittedStar; // To add to fsOutliers if it is a reference outlier.
90  if (measuredStar == nullptr) {
91  // it is a reference outlier
92  fittedStar = std::dynamic_pointer_cast<FittedStar>(chi2->star);
93  if (fittedStar->getMeasurementCount() == 0) {
94  LOGLS_WARN(_log, "FittedStar with no measuredStars found as an outlier: "
95  << *fittedStar << " chi2: " << chi2->chi2);
96  continue;
97  }
98  if (_nStarParams == 0) {
100  "RefStar is outlier but not removed when not fitting FittedStar-RefStar values: "
101  << *(fittedStar->getRefStar()) << " chi2: " << chi2->chi2);
102  continue;
103  }
104  // NOTE: Stars contribute twice to astrometry (x,y), but once to photometry (flux),
105  // NOTE: but we only need to mark one index here because both will be removed with that star.
106  indices.push_back(fittedStar->getIndexInMatrix());
107  LOGLS_TRACE(_log, "Removing refStar " << *(fittedStar->getRefStar()) << " chi2: " << chi2->chi2);
108  /* One might think it would be useful to account for PM
109  parameters here, but it is just useless */
110  } else {
111  // it is a measurement outlier
112  auto tempFittedStar = measuredStar->getFittedStar();
113  if (tempFittedStar->getMeasurementCount() == 1 && tempFittedStar->getRefStar() == nullptr) {
114  LOGLS_WARN(_log, "FittedStar with 1 measuredStar and no refStar found as an outlier: "
115  << *tempFittedStar);
116  continue;
117  }
118  getIndicesOfMeasuredStar(*measuredStar, indices);
119  LOGLS_TRACE(_log, "Removing measStar " << *measuredStar << " chi2: " << chi2->chi2);
120  }
121 
122  /* Find out if we already discarded a stronger outlier
123  constraining some parameter this one constrains as well. If
124  yes, we keep this one, because this stronger outlier could be
125  causing the large chi2 we have in hand. */
126  bool drop_it = true;
127  for (auto const &i : indices) {
128  if (affectedParams(i) != 0) {
129  drop_it = false;
130  }
131  }
132 
133  if (drop_it) // store the outlier in one of the lists:
134  {
135  if (measuredStar == nullptr) {
136  // reference term
137  fsOutliers.push_back(fittedStar);
138  } else {
139  // measurement term
140  msOutliers.push_back(measuredStar);
141  }
142  // mark the parameters as directly changed when we discard this chi2 term.
143  for (auto const &i : indices) {
144  affectedParams(i)++;
145  }
146  nOutliers++;
147  }
148  } // end loop on measurements/references
149  LOGLS_DEBUG(_log, "findOutliers: found " << msOutliers.size() << " meas outliers and "
150  << fsOutliers.size() << " ref outliers ");
151 
152  return nOutliers;
153 }
154 
155 namespace {
157 SparseMatrixD createHessian(std::size_t nParTot, TripletList const &tripletList) {
158  SparseMatrixD jacobian(nParTot, tripletList.getNextFreeIndex());
159  jacobian.setFromTriplets(tripletList.begin(), tripletList.end());
160  return jacobian * jacobian.transpose();
161 }
162 
164 void dumpMatrixAndGradient(SparseMatrixD const &matrix, Eigen::VectorXd const &grad,
165  std::string const &dumpFile, LOG_LOGGER _log) {
166  std::string ext = ".txt";
167  Eigen::MatrixXd matrixDense(matrix);
168  std::string dumpMatrixPath = dumpFile + "-mat" + ext;
169  std::ofstream matfile(dumpMatrixPath);
170  matfile << matrixDense << std::endl;
171  std::string dumpGradPath = dumpFile + "-grad" + ext;
172  std::ofstream gradfile(dumpGradPath);
173  gradfile << grad << std::endl;
174  LOGLS_INFO(_log, "Dumped Hessian, gradient to: '" << dumpMatrixPath << "', '" << dumpGradPath << "'");
175 }
176 } // namespace
177 
178 MinimizeResult FitterBase::minimize(std::string const &whatToFit, double nSigmaCut,
179  double sigmaRelativeTolerance, bool doRankUpdate, bool const doLineSearch,
180  std::string const &dumpMatrixFile) {
181  assignIndices(whatToFit);
182 
184 
185  // For the initial vector size, use all measured stars + all fitted stars, which should give the
186  // maximum possible number of triplets.
187  std::size_t nTrip = (_lastNTrip)
188  ? _lastNTrip
189  : _associations->getMaxMeasuredStars() + _associations->fittedStarList.size();
190  TripletList tripletList(nTrip);
191  Eigen::VectorXd grad(_nTotal);
192  grad.setZero();
193  double scale = 1.0;
194 
195  // Fill the triplets
196  leastSquareDerivatives(tripletList, grad);
197  _lastNTrip = tripletList.size();
198 
199  LOGLS_DEBUG(_log, "End of triplet filling, ntrip = " << tripletList.size());
200 
201  SparseMatrixD hessian = createHessian(_nTotal, tripletList);
202  tripletList.clear(); // we don't need it any more after we have the hessian.
203 
204  LOGLS_DEBUG(_log, "Starting factorization, hessian: dim="
205  << hessian.rows() << " non-zeros=" << hessian.nonZeros()
206  << " filling-frac = " << hessian.nonZeros() / std::pow(hessian.rows(), 2));
207 
208  if (dumpMatrixFile != "") {
209  if (hessian.rows() * hessian.cols() > 2e8) {
210  LOGLS_WARN(_log, "Hessian matrix is too big to dump to file, with rows, columns: "
211  << hessian.rows() << ", " << hessian.cols());
212  } else {
213  dumpMatrixAndGradient(hessian, grad, dumpMatrixFile, _log);
214  }
215  }
216 
218  if (chol.info() != Eigen::Success) {
219  LOGLS_ERROR(_log, "minimize: factorization failed ");
220  return MinimizeResult::Failed;
221  }
222 
223  std::size_t totalMeasOutliers = 0;
224  std::size_t totalRefOutliers = 0;
225  double oldChi2 = computeChi2().chi2;
226  double oldSigmaCut = 0;
227  double sigmaCut;
228 
229  while (true) {
230  Eigen::VectorXd delta = chol.solve(grad);
231  if (doLineSearch) {
232  scale = _lineSearch(delta);
233  }
234  offsetParams(scale * delta);
235  Chi2Statistic currentChi2(computeChi2());
236  LOGLS_DEBUG(_log, currentChi2);
237  if (!isfinite(currentChi2.chi2)) {
238  LOGL_ERROR(_log, "chi2 is not finite. Aborting outlier rejection.");
239  returnCode = MinimizeResult::NonFinite;
240  break;
241  }
242  if (currentChi2.chi2 > oldChi2 && totalMeasOutliers + totalRefOutliers != 0) {
243  LOGL_WARN(_log, "chi2 went up, skipping outlier rejection loop");
244  returnCode = MinimizeResult::Chi2Increased;
245  break;
246  }
247  oldChi2 = currentChi2.chi2;
248 
249  if (nSigmaCut == 0) break; // no rejection step to perform
250  MeasuredStarList msOutliers;
251  FittedStarList fsOutliers;
252  // keep nOutliers so we don't have to sum msOutliers.size()+fsOutliers.size() twice below.
253  std::size_t nOutliers = findOutliers(nSigmaCut, msOutliers, fsOutliers, sigmaCut);
254  double relChange = 1 - sigmaCut / oldSigmaCut;
255  LOGLS_DEBUG(_log, "findOutliers chi2 cut level: " << sigmaCut << ", relative change: " << relChange);
256  // If sigmaRelativeTolerance is set and at least one iteration has been done, break loop when the
257  // fractional change in sigmaCut levels is less than the sigmaRelativeTolerance parameter.
258  if ((sigmaRelativeTolerance > 0) && (oldSigmaCut > 0 && relChange < sigmaRelativeTolerance)) {
259  LOGLS_INFO(_log, "Iterations stopped with chi2 cut at " << sigmaCut << " and relative change of "
260  << relChange);
261  break;
262  }
263  totalMeasOutliers += msOutliers.size();
264  totalRefOutliers += fsOutliers.size();
265  oldSigmaCut = sigmaCut;
266  if (nOutliers == 0) break;
267  TripletList outlierTriplets(nOutliers);
268  grad.setZero(); // recycle the gradient
269  // compute the contributions of outliers to derivatives
270  outliersContributions(msOutliers, fsOutliers, outlierTriplets, grad);
271  // Remove significant outliers
272  removeMeasOutliers(msOutliers);
273  removeRefOutliers(fsOutliers);
274  if (doRankUpdate) {
275  // convert triplet list to eigen internal format
276  SparseMatrixD H(_nTotal, outlierTriplets.getNextFreeIndex());
277  H.setFromTriplets(outlierTriplets.begin(), outlierTriplets.end());
278  chol.update(H, false /* means downdate */);
279  // The contribution of outliers to the gradient is the opposite
280  // of the contribution of all other terms, because they add up to 0
281  grad *= -1;
282  } else {
283  // don't reuse tripletList because we want a new nextFreeIndex.
284  TripletList nextTripletList(_lastNTrip);
285  grad.setZero();
286  // Rebuild the matrix and gradient
287  leastSquareDerivatives(nextTripletList, grad);
288  _lastNTrip = nextTripletList.size();
289  LOGLS_DEBUG(_log, "Triplets recomputed, ntrip = " << nextTripletList.size());
290 
291  hessian = createHessian(_nTotal, nextTripletList);
292  nextTripletList.clear(); // we don't need it any more after we have the hessian.
293 
295  "Restarting factorization, hessian: dim="
296  << hessian.rows() << " non-zeros=" << hessian.nonZeros()
297  << " filling-frac = " << hessian.nonZeros() / std::pow(hessian.rows(), 2));
298  chol.compute(hessian);
299  if (chol.info() != Eigen::Success) {
300  LOGLS_ERROR(_log, "minimize: factorization failed ");
301  return MinimizeResult::Failed;
302  }
303  }
304  }
305 
306  if (totalMeasOutliers + totalRefOutliers > 0) {
307  _associations->cleanFittedStars();
308  }
309 
310  // only print the outlier summary if outlier rejection was turned on.
311  if (nSigmaCut != 0) {
312  LOGLS_INFO(_log, "Number of outliers (Measured + Reference = Total): "
313  << totalMeasOutliers << " + " << totalRefOutliers << " = "
314  << totalMeasOutliers + totalRefOutliers);
315  }
316  return returnCode;
317 }
318 
320  TripletList &tripletList, Eigen::VectorXd &grad) {
321  for (auto &outlier : msOutliers) {
322  MeasuredStarList tmp;
323  tmp.push_back(outlier);
324  const CcdImage &ccdImage = outlier->getCcdImage();
325  leastSquareDerivativesMeasurement(ccdImage, tripletList, grad, &tmp);
326  }
327  leastSquareDerivativesReference(fsOutliers, tripletList, grad);
328 }
329 
331  for (auto &measuredStar : outliers) {
332  auto fittedStar = measuredStar->getFittedStar();
333  measuredStar->setValid(false);
334  fittedStar->getMeasurementCount()--; // could be put in setValid
335  }
336 }
337 
339  for (auto &fittedStar : outliers) {
340  fittedStar->setRefStar(nullptr);
341  }
342 }
343 
344 void FitterBase::leastSquareDerivatives(TripletList &tripletList, Eigen::VectorXd &grad) const {
345  auto ccdImageList = _associations->getCcdImageList();
346  for (auto const &ccdImage : ccdImageList) {
347  leastSquareDerivativesMeasurement(*ccdImage, tripletList, grad);
348  }
349  leastSquareDerivativesReference(_associations->fittedStarList, tripletList, grad);
350 }
351 
352 void FitterBase::saveChi2Contributions(std::string const &baseName) const {
353  std::string replaceStr = "{type}";
354  auto pos = baseName.find(replaceStr);
355  std::string measFilename(baseName);
356  measFilename.replace(pos, replaceStr.size(), "-meas.csv");
357  std::string refFilename(baseName);
358  refFilename.replace(pos, replaceStr.size(), "-ref.csv");
359  saveChi2MeasContributions(measFilename);
360  saveChi2RefContributions(refFilename);
361 }
362 
363 double FitterBase::_lineSearch(Eigen::VectorXd const &delta) {
364  auto func = [this, &delta](double scale) {
365  auto offset = scale * delta;
366  offsetParams(offset);
367  auto chi2 = computeChi2();
368  // reset the system to where it was before offsetting.
369  offsetParams(-offset);
370  return chi2.chi2;
371  };
372  // The maximum theoretical precision is half the number of bits in the mantissa (see boost docs).
373  auto bits = std::numeric_limits<double>::digits / 2;
374  auto result = boost::math::tools::brent_find_minima(func, -1.0, 2.0, bits);
375  LOGLS_DEBUG(_log, "Line search scale factor: " << result.first);
376  return result.first;
377 }
378 
379 } // namespace jointcal
380 } // namespace lsst
py::object result
Definition: _schema.cc:429
Eigen::SparseMatrix< double, 0, Eigen::Index > SparseMatrixD
Definition: Eigenstuff.h:35
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 LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition: Log.h:547
#define LOGLS_INFO(logger, message)
Log a info-level message using an iostream-based interface.
Definition: Log.h:639
#define LOG_LOGGER
Definition: Log.h:714
#define LOGLS_ERROR(logger, message)
Log a error-level message using an iostream-based interface.
Definition: Log.h:679
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
Definition: Log.h:563
#define LOGLS_DEBUG(logger, message)
Log a debug-level message using an iostream-based interface.
Definition: Log.h:619
#define LOGLS_TRACE(logger, message)
Log a trace-level message using an iostream-based interface.
Definition: Log.h:599
T begin(T... args)
void update(SparseMatrixD const &H, bool UpOrDown)
Definition: Eigenstuff.h:68
Handler of an actual image from a single CCD.
Definition: CcdImage.h:64
Structure to accumulate the chi2 contributions per each star (to help find outliers).
Definition: Chi2.h:100
std::pair< double, double > computeAverageAndSigma()
Compute the average and std-deviation of these chisq values.
Definition: Chi2.cc:33
Simple structure to accumulate chi2 and ndof.
Definition: Chi2.h:52
A list of FittedStar s. Such a list is typically constructed by Associations.
Definition: FittedStar.h:116
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
void removeRefOutliers(FittedStarList &outliers)
Remove refStar outliers from the fit. No Refit done.
Definition: FitterBase.cc:338
virtual void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const =0
Set the indices of a measured star from the full matrix, for outlier removal.
Chi2Statistic computeChi2() const
Returns the chi2 for the current state.
Definition: FitterBase.cc:42
virtual void saveChi2MeasContributions(std::string const &filename) const =0
Save a CSV file containing residuals of measurement terms.
virtual void leastSquareDerivativesReference(FittedStarList const &fittedStarList, TripletList &tripletList, Eigen::VectorXd &grad) const =0
Compute the derivatives of the reference terms.
virtual void saveChi2Contributions(std::string const &baseName) const
Save the full chi2 term per star that was used in the minimization, for debugging.
Definition: FitterBase.cc:352
virtual void assignIndices(std::string const &whatToFit)=0
Set parameters to fit and assign indices in the big matrix.
virtual void offsetParams(Eigen::VectorXd const &delta)=0
Offset the parameters by the requested quantities.
std::shared_ptr< Associations > _associations
Definition: FitterBase.h:163
virtual void accumulateStatRefStars(Chi2Accumulator &accum) const =0
Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.
virtual void saveChi2RefContributions(std::string const &filename) const =0
Save a CSV file containing residuals of reference terms.
virtual void leastSquareDerivativesMeasurement(CcdImage const &ccdImage, TripletList &tripletList, Eigen::VectorXd &grad, MeasuredStarList const *measuredStarList=nullptr) const =0
Compute the derivatives of the measured stars and model for one CcdImage.
void outliersContributions(MeasuredStarList &msOutliers, FittedStarList &fsOutliers, TripletList &tripletList, Eigen::VectorXd &grad)
Contributions to derivatives from (presumably) outlier terms.
Definition: FitterBase.cc:319
std::size_t findOutliers(double nSigmaCut, MeasuredStarList &msOutliers, FittedStarList &fsOutliers, double &cut) const
Find Measurements and references contributing more than a cut, computed as.
Definition: FitterBase.cc:52
virtual void accumulateStatImageList(CcdImageList const &ccdImageList, Chi2Accumulator &accum) const =0
Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for measurements.
void removeMeasOutliers(MeasuredStarList &outliers)
Remove measuredStar outliers from the fit. No Refit done.
Definition: FitterBase.cc:330
MinimizeResult minimize(std::string const &whatToFit, double const nSigmaCut=0, double sigmaRelativeTolerance=0, bool const doRankUpdate=true, bool const doLineSearch=false, std::string const &dumpMatrixFile="")
Does a 1 step minimization, assuming a linear model.
Definition: FitterBase.cc:178
A list of MeasuredStar. They are usually filled in Associations::createCcdImage.
Definition: MeasuredStar.h:151
Eigen::Index getNextFreeIndex() const
Definition: Tripletlist.h:47
T clear(T... args)
T end(T... args)
T endl(T... args)
T find(T... args)
T isfinite(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
MinimizeResult
Return value of minimize()
Definition: FitterBase.h:40
A base class for image defects.
T pow(T... args)
T push_back(T... args)
T rbegin(T... args)
T rend(T... args)
T replace(T... args)
T reserve(T... args)
T size(T... args)
T sort(T... args)