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