LSST Applications g063fba187b+cac8b7c890,g0f08755f38+6aee506743,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+b4475c5878,g1dcb35cd9c+8f9bc1652e,g20f6ffc8e0+6aee506743,g217e2c1bcf+73dee94bd0,g28da252d5a+1f19c529b9,g2bbee38e9b+3f2625acfc,g2bc492864f+3f2625acfc,g3156d2b45e+6e55a43351,g32e5bea42b+1bb94961c2,g347aa1857d+3f2625acfc,g35bb328faa+a8ce1bb630,g3a166c0a6a+3f2625acfc,g3e281a1b8c+c5dd892a6c,g3e8969e208+a8ce1bb630,g414038480c+5927e1bc1e,g41af890bb2+8a9e676b2a,g7af13505b9+809c143d88,g80478fca09+6ef8b1810f,g82479be7b0+f568feb641,g858d7b2824+6aee506743,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+2903d499ea,gb58c049af0+d64f4d3760,gc28159a63d+3f2625acfc,gcab2d0539d+b12535109e,gcf0d15dbbd+46a3f46ba9,gda6a2b7d83+46a3f46ba9,gdaeeff99f8+1711a396fd,ge79ae78c31+3f2625acfc,gef2f8181fd+0a71e47438,gf0baf85859+c1f95f4921,gfa517265be+6aee506743,gfa999e8aa5+17cd334064,w.2024.51
LSST Data Management Base Package
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
lsst::jointcal::FitterBase Class Referenceabstract

Base class for fitters. More...

#include <FitterBase.h>

Inheritance diagram for lsst::jointcal::FitterBase:
lsst::jointcal::AstrometryFit lsst::jointcal::PhotometryFit

Public Member Functions

 FitterBase (std::shared_ptr< Associations > associations)
 
 FitterBase (FitterBase const &)=delete
 No copy or move: there is only ever one fitter of a given type.
 
 FitterBase (FitterBase &&)=delete
 
FitterBaseoperator= (FitterBase const &)=delete
 
FitterBaseoperator= (FitterBase &&)=delete
 
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.
 
Chi2Statistic computeChi2 () const
 Returns the chi2 for the current state.
 
void leastSquareDerivatives (TripletList &tripletList, Eigen::VectorXd &grad) const
 Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting.
 
virtual void offsetParams (Eigen::VectorXd const &delta)=0
 Offset the parameters by the requested quantities.
 
virtual void assignIndices (std::string const &whatToFit)=0
 Set parameters to fit and assign indices in the big matrix.
 
virtual void saveChi2Contributions (std::string const &baseName) const
 Save the full chi2 term per star that was used in the minimization, for debugging.
 
virtual ~FitterBase ()=default
 

Protected Member Functions

virtual void saveChi2MeasContributions (std::string const &filename) const =0
 Save a CSV file containing residuals of measurement terms.
 
virtual void saveChi2RefContributions (std::string const &filename) const =0
 Save a CSV file containing residuals of reference 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.
 
void outliersContributions (MeasuredStarList &msOutliers, FittedStarList &fsOutliers, TripletList &tripletList, Eigen::VectorXd &grad)
 Contributions to derivatives from (presumably) outlier terms.
 
void removeMeasOutliers (MeasuredStarList &outliers)
 Remove measuredStar outliers from the fit. No Refit done.
 
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.
 
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.
 
virtual void accumulateStatRefStars (Chi2Accumulator &accum) const =0
 Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.
 
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.
 
virtual void leastSquareDerivativesReference (FittedStarList const &fittedStarList, TripletList &tripletList, Eigen::VectorXd &grad) const =0
 Compute the derivatives of the reference terms.
 

Protected Attributes

std::shared_ptr< Associations_associations
 
std::string _whatToFit
 
Eigen::Index _lastNTrip
 
Eigen::Index _nTotal
 
Eigen::Index _nModelParams
 
Eigen::Index _nStarParams
 
LOG_LOGGER _log
 

Detailed Description

Base class for fitters.

Implements minimize and findOutliers. Chi2, residual, derivative, etc. calculations must be implemented in the child class via the virtual methods.

Definition at line 55 of file FitterBase.h.

Constructor & Destructor Documentation

◆ FitterBase() [1/3]

lsst::jointcal::FitterBase::FitterBase ( std::shared_ptr< Associations > associations)
inlineexplicit

Definition at line 57 of file FitterBase.h.

58 : _associations(std::move(associations)),
59 _whatToFit(""),
60 _lastNTrip(0),
61 _nTotal(0),
63 _nStarParams(0) {}
std::shared_ptr< Associations > _associations
Definition FitterBase.h:165
T move(T... args)

◆ FitterBase() [2/3]

lsst::jointcal::FitterBase::FitterBase ( FitterBase const & )
delete

No copy or move: there is only ever one fitter of a given type.

◆ FitterBase() [3/3]

lsst::jointcal::FitterBase::FitterBase ( FitterBase && )
delete

◆ ~FitterBase()

virtual lsst::jointcal::FitterBase::~FitterBase ( )
virtualdefault

Member Function Documentation

◆ accumulateStatImageList()

virtual void lsst::jointcal::FitterBase::accumulateStatImageList ( CcdImageList const & ccdImageList,
Chi2Accumulator & accum ) const
protectedpure virtual

Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for measurements.

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

◆ accumulateStatRefStars()

virtual void lsst::jointcal::FitterBase::accumulateStatRefStars ( Chi2Accumulator & accum) const
protectedpure virtual

Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

◆ assignIndices()

virtual void lsst::jointcal::FitterBase::assignIndices ( std::string const & whatToFit)
pure virtual

Set parameters to fit and assign indices in the big matrix.

Parameters
[in]whatToFitSee child class documentation for valid string values.

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

◆ computeChi2()

Chi2Statistic lsst::jointcal::FitterBase::computeChi2 ( ) const

Returns the chi2 for the current state.

Definition at line 43 of file FitterBase.cc.

43 {
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}
virtual void accumulateStatRefStars(Chi2Accumulator &accum) const =0
Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars.
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.

◆ findOutliers()

std::size_t lsst::jointcal::FitterBase::findOutliers ( double nSigmaCut,
MeasuredStarList & msOutliers,
FittedStarList & fsOutliers,
double & cut ) const
protected

Find Measurements and references contributing more than a cut, computed as.

\[ <chi2> + nSigmaCut + rms(chi2). \]

The outliers are NOT removed, and no refit is done.

After returning from here, there are still measurements that contribute above the cut, but their contribution should be evaluated after a refit before discarding them.

Parameters
[in]nSigmaCutNumber of sigma to select on.
[out]msOutlierslist of MeasuredStar outliers to populate
[out]fsOutlierslist of FittedStar outliers to populate
[out]cutvalue of chi2 that defines which objects are outliers
Returns
Total number of outliers that were removed.

Definition at line 53 of file FitterBase.cc.

54 {
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}
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition Log.h:659
#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
virtual void getIndicesOfMeasuredStar(MeasuredStar const &measuredStar, IndexVector &indices) const =0
Set the indices of a measured star from the full matrix, for outlier removal.
T push_back(T... args)
T sort(T... args)

◆ getIndicesOfMeasuredStar()

virtual void lsst::jointcal::FitterBase::getIndicesOfMeasuredStar ( MeasuredStar const & measuredStar,
IndexVector & indices ) const
protectedpure virtual

Set the indices of a measured star from the full matrix, for outlier removal.

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

◆ leastSquareDerivatives()

void lsst::jointcal::FitterBase::leastSquareDerivatives ( TripletList & tripletList,
Eigen::VectorXd & grad ) const

Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting.

The Jacobian is given as triplets in a sparse matrix, the gradient as a dense vector. The parameters which vary, and their indices, are to be set using assignIndices.

Parameters
tripletListtripletList of (row,col,value) representing the Jacobian of the chi2.
gradThe gradient of the chi2.

Definition at line 347 of file FitterBase.cc.

347 {
348 auto ccdImageList = _associations->getCcdImageList();
349 for (auto const &ccdImage : ccdImageList) {
350 leastSquareDerivativesMeasurement(*ccdImage, tripletList, grad);
351 }
352 leastSquareDerivativesReference(_associations->fittedStarList, tripletList, grad);
353}
virtual void leastSquareDerivativesReference(FittedStarList const &fittedStarList, TripletList &tripletList, Eigen::VectorXd &grad) const =0
Compute the derivatives of the 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.

◆ leastSquareDerivativesMeasurement()

virtual void lsst::jointcal::FitterBase::leastSquareDerivativesMeasurement ( CcdImage const & ccdImage,
TripletList & tripletList,
Eigen::VectorXd & grad,
MeasuredStarList const * measuredStarList = nullptr ) const
protectedpure virtual

Compute the derivatives of the measured stars and model for one CcdImage.

The last argument will process a sub-list for outlier removal.

Implemented in lsst::jointcal::PhotometryFit, and lsst::jointcal::AstrometryFit.

◆ leastSquareDerivativesReference()

virtual void lsst::jointcal::FitterBase::leastSquareDerivativesReference ( FittedStarList const & fittedStarList,
TripletList & tripletList,
Eigen::VectorXd & grad ) const
protectedpure virtual

Compute the derivatives of the reference terms.

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

◆ minimize()

MinimizeResult lsst::jointcal::FitterBase::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.

This is a complete Newton Raphson step. Compute first and second derivatives, solve for the step and apply it, with an optional line search.

It calls assignIndices, leastSquareDerivatives, solves the linear system and calls offsetParams, then removes outliers in a loop if requested. Relies on sparse linear algebra via Eigen's CholmodSupport package.

Parameters
[in]whatToFitSee child method assignIndices for valid string values.
[in]nSigmaCutHow many sigma to reject outliers at. Outlier rejection ignored for nSigmaCut=0.
[in]sigmaRelativeTolerancePercentage change in the chi2 cut for outliers tolerated for termination. If value is zero, minimization iterations will continue until there are no outliers.
[in]doRankUpdateUse CholmodSimplicialLDLT2.update() to do a fast rank update after outlier removal; otherwise do a slower full recomputation of the matrix. Only matters if nSigmaCut != 0.
[in]doLineSearchUse boost's brent_find_minima to perform a line search after the gradient solution is found, and apply the scale factor to the computed offsets. The line search is done in the domain [-1, 2], but if the scale factor is far from 1.0, then the problem is likely in a significantly non-linear regime.
[in]dumpMatrixFileWrite the pre-fit Hessian matrix and gradient to the files with "-mat.txt" and "-grad.txt". Be aware, this requires a large increase in memory usage to create a dense matrix before writing it; the output file may be large. Writing the matrix can be helpful for debugging bad fits. Read it and compute the real eigenvalues (recall that the Hessian is symmetric by construction) with numpy:
hessian = np.matrix(np.loadtxt("dumpMatrixFile-mat.txt"))
values, vectors = np.linalg.eigh(hessian)
Returns
Return code describing success/failure of fit.
Note
When fitting one parameter set by itself (e.g. "Model"), the system is purely linear (assuming there are no cross-terms in the derivatives, e.g. the SimpleAstrometryModel), which should result in the optimal chi2 after a single step. This can be used to debug the fitter by fitting that parameter set twice in a row: the second run with the same "whatToFit" will produce no change in the fitted parameters, if the calculations and indices are defined correctly.

Definition at line 179 of file FitterBase.cc.

181 {
182 assignIndices(whatToFit);
183
184 MinimizeResult returnCode = MinimizeResult::Converged;
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}
Eigen::SparseMatrix< double, 0, Eigen::Index > SparseMatrixD
Definition Eigenstuff.h:35
#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 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
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.
Chi2Statistic computeChi2() const
Returns the chi2 for the current state.
Definition FitterBase.cc:43
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.
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
void removeMeasOutliers(MeasuredStarList &outliers)
Remove measuredStar outliers from the fit. No Refit done.
T isfinite(T... args)
scale(algorithm, min, max=None, frame=None)
Definition ds9.py:108
T pow(T... args)

◆ offsetParams()

virtual void lsst::jointcal::FitterBase::offsetParams ( Eigen::VectorXd const & delta)
pure virtual

Offset the parameters by the requested quantities.

The used parameter layout is the one from the last call to assignIndices or minimize(). There is no easy way to check that the current setting of whatToFit and the provided Delta vector are compatible: we can only test the size.

Parameters
[in]deltavector of offsets to apply

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

◆ operator=() [1/2]

FitterBase & lsst::jointcal::FitterBase::operator= ( FitterBase && )
delete

◆ operator=() [2/2]

FitterBase & lsst::jointcal::FitterBase::operator= ( FitterBase const & )
delete

◆ outliersContributions()

void lsst::jointcal::FitterBase::outliersContributions ( MeasuredStarList & msOutliers,
FittedStarList & fsOutliers,
TripletList & tripletList,
Eigen::VectorXd & grad )
protected

Contributions to derivatives from (presumably) outlier terms.

No discarding done.

Definition at line 322 of file FitterBase.cc.

323 {
324 for (auto &outlier : msOutliers) {
325 MeasuredStarList tmp;
326 tmp.push_back(outlier);
327 const CcdImage &ccdImage = outlier->getCcdImage();
328 leastSquareDerivativesMeasurement(ccdImage, tripletList, grad, &tmp);
329 }
330 leastSquareDerivativesReference(fsOutliers, tripletList, grad);
331}

◆ removeMeasOutliers()

void lsst::jointcal::FitterBase::removeMeasOutliers ( MeasuredStarList & outliers)
protected

Remove measuredStar outliers from the fit. No Refit done.

Definition at line 333 of file FitterBase.cc.

333 {
334 for (auto &measuredStar : outliers) {
335 auto fittedStar = measuredStar->getFittedStar();
336 measuredStar->setValid(false);
337 fittedStar->getMeasurementCount()--; // could be put in setValid
338 }
339}

◆ removeRefOutliers()

void lsst::jointcal::FitterBase::removeRefOutliers ( FittedStarList & outliers)
protected

Remove refStar outliers from the fit. No Refit done.

Definition at line 341 of file FitterBase.cc.

341 {
342 for (auto &fittedStar : outliers) {
343 fittedStar->setRefStar(nullptr);
344 }
345}

◆ saveChi2Contributions()

void lsst::jointcal::FitterBase::saveChi2Contributions ( std::string const & baseName) const
virtual

Save the full chi2 term per star that was used in the minimization, for debugging.

Saves results to text files "baseName-meas.csv" and "baseName-ref.csv" for the MeasuredStar and RefStar contributions, respectively. This method is mostly useful for debugging: we will probably want to create a better persistence system for jointcal's internal representations in the future (see DM-12446).

Definition at line 355 of file FitterBase.cc.

355 {
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}
virtual void saveChi2MeasContributions(std::string const &filename) const =0
Save a CSV file containing residuals of measurement terms.
virtual void saveChi2RefContributions(std::string const &filename) const =0
Save a CSV file containing residuals of reference terms.
T find(T... args)
T size(T... args)

◆ saveChi2MeasContributions()

virtual void lsst::jointcal::FitterBase::saveChi2MeasContributions ( std::string const & filename) const
protectedpure virtual

Save a CSV file containing residuals of measurement terms.

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

◆ saveChi2RefContributions()

virtual void lsst::jointcal::FitterBase::saveChi2RefContributions ( std::string const & filename) const
protectedpure virtual

Save a CSV file containing residuals of reference terms.

Implemented in lsst::jointcal::AstrometryFit, and lsst::jointcal::PhotometryFit.

Member Data Documentation

◆ _associations

std::shared_ptr<Associations> lsst::jointcal::FitterBase::_associations
protected

Definition at line 165 of file FitterBase.h.

◆ _lastNTrip

Eigen::Index lsst::jointcal::FitterBase::_lastNTrip
protected

Definition at line 168 of file FitterBase.h.

◆ _log

LOG_LOGGER lsst::jointcal::FitterBase::_log
protected

Definition at line 174 of file FitterBase.h.

◆ _nModelParams

Eigen::Index lsst::jointcal::FitterBase::_nModelParams
protected

Definition at line 170 of file FitterBase.h.

◆ _nStarParams

Eigen::Index lsst::jointcal::FitterBase::_nStarParams
protected

Definition at line 171 of file FitterBase.h.

◆ _nTotal

Eigen::Index lsst::jointcal::FitterBase::_nTotal
protected

Definition at line 169 of file FitterBase.h.

◆ _whatToFit

std::string lsst::jointcal::FitterBase::_whatToFit
protected

Definition at line 166 of file FitterBase.h.


The documentation for this class was generated from the following files: