LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
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. More...
 
 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. More...
 
Chi2Statistic computeChi2 () const
 Returns the chi2 for the current state. More...
 
void leastSquareDerivatives (TripletList &tripletList, Eigen::VectorXd &grad) const
 Evaluates the chI^2 derivatives (Jacobian and gradient) for the current whatToFit setting. More...
 
virtual void offsetParams (Eigen::VectorXd const &delta)=0
 Offset the parameters by the requested quantities. More...
 
virtual void assignIndices (std::string const &whatToFit)=0
 Set parameters to fit and assign indices in the big matrix. More...
 
virtual void saveChi2Contributions (std::string const &baseName) const
 Save the full chi2 term per star that was used in the minimization, for debugging. More...
 
virtual ~FitterBase ()=default
 

Protected Member Functions

virtual void saveChi2MeasContributions (std::string const &filename) const =0
 Save a CSV file containing residuals of measurement terms. More...
 
virtual void saveChi2RefContributions (std::string const &filename) const =0
 Save a CSV file containing residuals of reference terms. More...
 
std::size_t findOutliers (double nSigmaCut, MeasuredStarList &msOutliers, FittedStarList &fsOutliers, double &cut) const
 Find Measurements and references contributing more than a cut, computed as. More...
 
void outliersContributions (MeasuredStarList &msOutliers, FittedStarList &fsOutliers, TripletList &tripletList, Eigen::VectorXd &grad)
 Contributions to derivatives from (presumably) outlier terms. More...
 
void removeMeasOutliers (MeasuredStarList &outliers)
 Remove measuredStar outliers from the fit. No Refit done. More...
 
void removeRefOutliers (FittedStarList &outliers)
 Remove refStar outliers from the fit. No Refit done. More...
 
virtual void getIndicesOfMeasuredStar (MeasuredStar const &measuredStar, IndexVector &indices) const =0
 Set the indices of a measured star from the full matrix, for outlier removal. More...
 
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. More...
 
virtual void accumulateStatRefStars (Chi2Accumulator &accum) const =0
 Compute the chi2 (per star or total, depending on which Chi2Accumulator is used) for RefStars. More...
 
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. More...
 
virtual void leastSquareDerivativesReference (FittedStarList const &fittedStarList, TripletList &tripletList, Eigen::VectorXd &grad) const =0
 Compute the derivatives of the reference terms. More...
 

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) {}
Eigen::Index _nModelParams
Definition: FitterBase.h:170
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.

◆ 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.

◆ 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 42 of file FitterBase.cc.

42 {
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}
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 52 of file FitterBase.cc.

53 {
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}
#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.

◆ 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 346 of file FitterBase.cc.

346 {
347 auto ccdImageList = _associations->getCcdImageList();
348 for (auto const &ccdImage : ccdImageList) {
349 leastSquareDerivativesMeasurement(*ccdImage, tripletList, grad);
350 }
351 leastSquareDerivativesReference(_associations->fittedStarList, tripletList, grad);
352}
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.

◆ 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.

◆ 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 178 of file FitterBase.cc.

180 {
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 ");
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 = 0.;
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");
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 = 0.;
255 if(oldSigmaCut!=0.) relChange = (1 - sigmaCut / oldSigmaCut);
256
257 LOGLS_DEBUG(_log, "findOutliers chi2 cut level: " << sigmaCut << ", relative change: " << relChange);
258 // If sigmaRelativeTolerance is set and at least one iteration has been done, break loop when the
259 // fractional change in sigmaCut levels is less than the sigmaRelativeTolerance parameter.
260 if ((sigmaRelativeTolerance > 0) && (oldSigmaCut > 0 && relChange < sigmaRelativeTolerance)) {
261 LOGLS_INFO(_log, "Iterations stopped with chi2 cut at " << sigmaCut << " and relative change of "
262 << relChange);
263 break;
264 }
265 totalMeasOutliers += msOutliers.size();
266 totalRefOutliers += fsOutliers.size();
267 oldSigmaCut = sigmaCut;
268 if (nOutliers == 0) break;
269 TripletList outlierTriplets(nOutliers);
270 grad.setZero(); // recycle the gradient
271 // compute the contributions of outliers to derivatives
272 outliersContributions(msOutliers, fsOutliers, outlierTriplets, grad);
273 // Remove significant outliers
274 removeMeasOutliers(msOutliers);
275 removeRefOutliers(fsOutliers);
276 if (doRankUpdate) {
277 // convert triplet list to eigen internal format
278 SparseMatrixD H(_nTotal, outlierTriplets.getNextFreeIndex());
279 H.setFromTriplets(outlierTriplets.begin(), outlierTriplets.end());
280 chol.update(H, false /* means downdate */);
281 // The contribution of outliers to the gradient is the opposite
282 // of the contribution of all other terms, because they add up to 0
283 grad *= -1;
284 } else {
285 // don't reuse tripletList because we want a new nextFreeIndex.
286 TripletList nextTripletList(_lastNTrip);
287 grad.setZero();
288 // Rebuild the matrix and gradient
289 leastSquareDerivatives(nextTripletList, grad);
290 _lastNTrip = nextTripletList.size();
291 LOGLS_DEBUG(_log, "Triplets recomputed, ntrip = " << nextTripletList.size());
292
293 hessian = createHessian(_nTotal, nextTripletList);
294 nextTripletList.clear(); // we don't need it any more after we have the hessian.
295
297 "Restarting factorization, hessian: dim="
298 << hessian.rows() << " non-zeros=" << hessian.nonZeros()
299 << " filling-frac = " << hessian.nonZeros() / std::pow(hessian.rows(), 2));
300 chol.compute(hessian);
301 if (chol.info() != Eigen::Success) {
302 LOGLS_ERROR(_log, "minimize: factorization failed ");
304 }
305 }
306 }
307
308 if (totalMeasOutliers + totalRefOutliers > 0) {
309 _associations->cleanFittedStars();
310 }
311
312 // only print the outlier summary if outlier rejection was turned on.
313 if (nSigmaCut != 0) {
314 LOGLS_INFO(_log, "Number of outliers (Measured + Reference = Total): "
315 << totalMeasOutliers << " + " << totalRefOutliers << " = "
316 << totalMeasOutliers + totalRefOutliers);
317 }
318 return returnCode;
319}
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.
Definition: FitterBase.cc:346
void removeRefOutliers(FittedStarList &outliers)
Remove refStar outliers from the fit. No Refit done.
Definition: FitterBase.cc:340
Chi2Statistic computeChi2() const
Returns the chi2 for the current state.
Definition: FitterBase.cc:42
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.
Definition: FitterBase.cc:321
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
void removeMeasOutliers(MeasuredStarList &outliers)
Remove measuredStar outliers from the fit. No Refit done.
Definition: FitterBase.cc:332
T isfinite(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
MinimizeResult
Return value of minimize()
Definition: FitterBase.h:42
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 321 of file FitterBase.cc.

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

◆ removeMeasOutliers()

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

Remove measuredStar outliers from the fit. No Refit done.

Definition at line 332 of file FitterBase.cc.

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

◆ removeRefOutliers()

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

Remove refStar outliers from the fit. No Refit done.

Definition at line 340 of file FitterBase.cc.

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

◆ 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 354 of file FitterBase.cc.

354 {
355 std::string replaceStr = "{type}";
356 auto pos = baseName.find(replaceStr);
357 std::string measFilename(baseName);
358 measFilename.replace(pos, replaceStr.size(), "-meas.csv");
359 std::string refFilename(baseName);
360 refFilename.replace(pos, replaceStr.size(), "-ref.csv");
361 saveChi2MeasContributions(measFilename);
362 saveChi2RefContributions(refFilename);
363}
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: